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A  SCORE  FOR  CORRECT  DATA  ASSOCIATION  IN  MULTI-TARGET  TRACKING 


D.  L  Ahpach  and  R.  N.  Lobbia 
ORINCON  CORPORAT  ION 
La  Jolla.  California  9203? 


SUMMARY 

la  the  real- world  multi-target  tracking  problem,  there  exists  the  potsi- 
bOity  for  many  thinp  to  go  wrong.  Typical  probiem*  which  ariae 
include:  too  few  tracks  are  formed;  too  many  tracks  are  formed  (false 
tracks);  and  inaccurate  position,  course,  and  speed  estimates  are 
reported.  The  above  difficulties  ate  often  the  result  of  incorrect  alloca¬ 
tion  of  data  to  individual  tracks.  Algorithms,  while  estimating  the 
motion  of  a  given  target,  inadvertently  mix  in  clutter  and/or  measure¬ 
ments  from  another  target,  la  order  for  correct  allocation  of  data  to  a 
given  track  to  be  made,  one  must  have  an  effective  scoring  formula; 
that  is,  some  means  of  determining  how  likely  a  given  assignment  of 
data  is.  To  be  effective,  a  scoring  formula  must  produce  (on  the  aver¬ 
age)  a  better  score  for  correct  assignments  than  for  incorrect  assign¬ 
ments.  Information  useful  in  the  scoring  process  includes  a  priori  intel¬ 
ligence  data  (such  as  initial  target  locations),  models  of  target  motion, 
models  of  the  transmission  channel,  and  expected  moments  of  clutter 
for  the  sensor  gain  setting  being  used.  Basically,  the  score  is  derived 
from  the  residuals  which  come  out  of  the  processing  of  a  batch  of  data 
with  the  extended  Kalman  filter.  This  it  used  to  evaluate  the  likelihood 
of  potential  tracks.  Although  the  “likelihood''  has  an  intuitive  meaning, 
the  term  is  used  here  to  mean  the  probability  density  function  p(  X)  of 
the  track  X.  The  expected  cost  of  a  given  assignment  is  derived  with 
the  theory  of  extremals  being  used  to  obtain  the  expected  cost  of  add¬ 
ing  a  clutter  point  in  a  track.  The  resulting  expected  cost  is  then  shown 
to  behave  in  a  quantitative  fashion  and  this  can  be  visualized  from  a 
geometric  viewpoint. 

1.  INTRODUCTION 

In  the  last  few  years  a  number  of  approaches  to  the  problem  of 
tracking  multiple  targets  in  a  cluttered  environment  have  been  pub¬ 
lished.  Some  aspects  of  the  problem  that  have  been  considered  in 
some  subeets  of  these  publications  include  the  problem  of  false  alarms 
or  missing  measurements,  track  initialization,  multiple  measurement 
types  (MSI)  and  target  classification. 

Many  of  the  so-called  multi-target  trackers  described  in  the  open 
literature  really  deal  only  with  the  problem  of  one  target  in  clutter. 

This  hypothesis  of  only  one  target  can  greatly  restrict  the  viability  of 
a  multi-target  tracker  to  sort  out  confused  situations. 

Of  the  tnckers  that  have  been  proposed,  and  in  some  cases  imple¬ 
mented,  one  can  see  certain  similarities  and  differences  which  allow 
the  trackers  to  be  grouped  into  certain  classes.  The  grouping  and 
certain  applications  of  these  groups  has  led  to  the  following  thoughts. 

Perhaps  the  most  fundamental  aspect  of  a  tracker  is  how  it  handles 
and  interacts  with  the  data.  The  data  is  after  all  our  handle  on  the  real 
world  and  all  the  information  we  have  about  a  specific  tracking  realiza¬ 
tion  is  contained  in  the  data.  A  second  fundamental  aspect  of  the 
multi-target  tracking  in  clutter  problem  a  that  of  alternate  hypotheses. 
It  is  possible  to  derive  the  probability  distribution  of  tracks  and  clutter 
points  if  one  carefully  specifies  the  a  priori  probability  base,  i.e.,  the 
probabilistic  target  models,  probabilistic  measurement  models,  etc. 

In  very  simple  cases  where  one  can  get  the  optimal  solution,  one  finds 
this  consists  of  all  possible  configurations  of  the  data  into  the  sets.  In 
each  configuration  each  data  set  represents  a  possible  alternate  track 
or  the  set  of  clutter  points.  A  probability  measure  is  assigned  to  each 
possible  total  surveillance  region  picture.  This  globally  optimal 
approach  is  generally  not  reasonable  for  implementation  and  approxi¬ 
mate  or  suboptima)  approaches  must  be  considered. 

Several  approaches  focus  entirely  on  the  construction  of  the  sets  of 
data  (the  construction  of  feasible  tracks).  It  is  this  latter  philosophy 
that  is  being  addressed  in  this  paper. 

Conceptually,  using  maximum  likelihood  techniques,  various  com¬ 
binations  of  data  are  tried  and  then  “scured"  using  log-likelihood 
functions.  The  best  fit  to  thc-urget  model  gets  the  lowest  score  and 
this  is  considered  to  be  one  of  the  tracks  in  the  region  if  no  other  low 


score  track  competes  for  the  same  measurements.  If  two  or  more 
targets  compete  for  the  tame  measurements,  several  situetiom  can 
occur.  These  include  the  poesibilities  that  the  two  targets  are  lumped 
together,  one  target  is  rejected,  the  targets  get  mixed  with  track  points 
assigned  to  clutter  and  two  “bad"  tracks  reported,  or  she  care  that  all 
points  ire  assigned  to  clutter.  It  is  pooibtc,  though  perhaps  not  nor¬ 
mal,  to  find  situations  where  the  choice  of  the  beat  nonovertappmg 
feasible  track  docs  not  correspond  to  a  best  surveillance  picture.  This 
it  quite  easy  to  do  if  the  tracks  overlap  and  compete  for  the  aaaae 
measurements. 

Many  trackers  consider  alternate  hypotheses  at  far  as  atrigniag  mea¬ 
surements  to  a  track.  However,  once  a  decision  hat  been  made  that  a 
measurement  belongs  with  another  group  of  measurements,  this  deci¬ 
sion  is  not  re-examined.  Once  the  decision  has  been  made  to  areign  a 
piece  of  data  to  a  “track,”  that  decision  is  final.  This  is  done  horante 
the  system  usually  requites  “an  answer.”  Also,  there  is  always  new 
data  coming  into  the  tracker  allowing  new  hypothesis  tests.  In  addi¬ 
tion,  there  is  a  limited  amount  of  computer  response.  One  could  say 
that  the  hypothesis  testing  is  directed  to  make  a  decision  on  the 
proper  surveillance  picture  or  that  the  tncker  it  “decision  directed.” 

In  the  next  section,  we  will  set  how  effective  scoring  algorithms  are 
developed -ones  which  can  handle  the  alternate  scenarios  posed  above. 
Following  this,  in  Section  3,  a  refinement  to  this  scoring  algorithm  is 
proposed  and  it  is  seen  that  the  avenge  cost  incurred  for  assigning 
measurements  to  tracks  can  be  visualized  from  a  geometrical  view- 
point.  In  particular,  it  will  be  shown  quantitatively  tint  there  exists 
a  unique  number  or  points  (measurements)  in  a  given  track  that  yields 
an  average  minimum  cost.  Incorrectly  assigning  clutter  points  to  this 
tnck  and  wrongly  assigning  points  to  duller  will,  on  the  average, 
increase  the  cost  in  a  well-defined  manner. 


2.  SCORINC  ALGORITHMS 

In  order  for  a  correct  assignment  of  measurement  data  to  a  given 
track  to  be  made,  we  must  have  an  effective  scoring  formula,  i.c., 
some  means  of  determining  how  likely  a  given  assignment  of  data  is. 

To  be  effective,  a  scoring  formula  mutt  produce  (on  the  avenge)  a 
better  score  for  correct  assignments  than  for  incorrect  assignments. 
Informition  useful  in  the  scoring  process  includes  a  priori  intelligence 
data  (such  as  inithl  target  locations),  models  of  target  motion,  and 
expected  amounts  of  clutter  for  the  sensor  gain  setting  being  used. 
Basically,  the  score  Is  derived  from  the  residuals  which  come  out  of 
the  processing  of  a  batch  of  data  with  the  extended  Kalman  filter.  This 
is  used  to  evaluate  the  likelihood  of  potential  tracks.  Although  “likeli¬ 
hood"  has  t  useful  intuitive  meaning,  we  me  the  term  to  mean  the 
probability  density  function  p(  X)  of  the  track  X.  The  concepts  we  we 
are  well-known,  since  most  of  the  work  in  estimation  theory  pertains 
to  situations  where  all  the  observations  Z  *  {t|,  la,  •  ■  •,  are  due  to 
a  single  target.  An  obvious  example  is  the  stochastic  linear  system 


*k+l 

•  *kxk  *  ®kuk> 

k«0.  l,---,n. 

(1) 

*k 

■  ckxk*¥k- 

k»  1.  •••.», 

(2) 

with  states  { xk)  C  R*,  observations  {zk)  C  R*.  process  noise  {  uk)  C  R 
and  measurement  noise  (vj,  ]  c  Rv.  Ak.  B|j,  Ck  are  matrices  of  appropri¬ 
ate  dimension  that  may  vary  with  time.  The  initial  stale  xq  is  a  Gaussian 
random  vector  with  covariance  Pq,  independent  of  the  processes  (uk) 
and  {vk},  which  art  themselves  zrro mean  white  Cauttian  noise  with 
covariances  (Qfc]  ami  (Kk).  Under  these  axMimptmm,  the  weB 
known  Kalman  equations  provide  minimum  variance  unbiased  estimates 
{ fife)  of  the  states  bated  on  alt  past  data; 


Ucftarr  wry  ikw  fxunl  k  otldrd  lu a  pailial  luki,  il  dumb!  |m»*  a 


Aktk«Kk(xko|4k«|Ak£k). 
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*kcUl<ck-rl\Ck-M  *  Kk+l)’1- 

<43 

VkAl*BkQk“x. 

IS) 

(l-KkCk+1)?k,  k  ■  0, 1.  n . 

(6) 

When  nonlinear  measurements  are  involved,  a  simple  linearization 
process (the  extended  Kalman  filler)  is  used.  The  equations  remain 
exactly  the  same,  except  that  the  term 

xk*rck»lAklk  <7) 


coarse  test, 

B*klL  <  Ok  <'7' 

a  fine  test, 

«ltVi,*k<7k  «'«» 

and,  finally,  the  likelihood  test.  The  coarse  test  checks  the  mapnitude  of 
the  maximum  component  of  the  vector  Sk  «  Rz  agamst  6k.  and  is 
included  because  it  is  computationally  cheaper  to  perform  than  the  fou¬ 
lest.  The  constants  4k  and  yk  can  bo  chosen  so  that 


i»  the  first  equation  above  is  replaced  by 

*kel'hke|<Ak*k>»  (g) 

where  hk+ j(-) denotes  the  nonlinear  relabionship  between  the  measure¬ 
ments  and  state  vector,  and  the  Cj, -matrix  becomes 


Ck 


ahk(x> 

dx 


Mk-i 


(9) 


It  is  natural  to  compute  the  likelihood  function  p(  X )  for  the  track 
X  C  Z  based  on  the  Kalman  filter  state  estimates.  The  innovations 
sequence  is  an  integral  part  of  this  computation,  which  is  a  sequence  of 
the  measurement  residuals: 

4k+l  *  ik+i  -  Ck+|  Akik  ,  k  “  0, 1,  — ,  n.  (10) 

The  (negative)  log  likelihood  function  is  given  in  terms  of  {4  k}  by 
n  n 

c(X)  »  ndim(z)ln  2»  *  £  ln|Vkl  *  -L  £  4jtVit  4|c  (ID 

*  k-1  *  k-l 


The  covariance  matrix  for  the  measurement  residual,  Vk,  can  be  com¬ 
puted  directly: 

Vk  -  Ck+1?kCj.+,  ♦  Rk+1.  k«  1. 2,  n  .  02) 


{*k  l*itv»t  «k  «  ’k(  {*kl  H«klL  <  *k}-  CW) 

One  difficulty  with  this  approach  is  that  because  of  the  ‘'deterministic" 
terms  in  the  likelihood  function 

n 

n  dm(z)  In  2*  +1  £  |n  |Vk|  (20) 

*  k«l 

it  is  difficult  to  compare  some  of  the  tracks  of  different  length.  It  is 
also  difficult  to  assess  the  absolute  goodness  of  a  score.  Therefore,  an 
alternate  score  with  a  more  absolute  meaning  can  be  defined.  This  is 


described  in  the  next  section. 

3.  REFINED  SCORING 

For  track  i  at  stage  k  define  the  stagewisc  chi-squared  score 

Sic  ■  (Mc'HjTv-M-H)  •  <21  > 

This  has  the  fcalurcs-ir  all  measurements  have  been  assigned  to  the 
correct  targets  and  all  filter  parameters  chosen  corfectly-that 

EJSy  •  E(Sk)  »  2  (22) 

E  {(Sk  -  2)2  )  ■  4  123) 

°Sk  “  2  (24) 


Each  feasible  track  is  the  result  of  a  hypothesis  test  that  uses  the 
track  likelihood  function  p(  X )  fur  equivalently,  the  negative  log  likeli¬ 
hood  c(X))  determined  from  the  Kalman  filter.  Since  the  density 
function  of  the  alternative  hypothesis  (that  X  is  not  a  track)  is  unknown. 


the  decision  rule  is  simply 

-«X)  ■  In p(X  |  {Sk )(ja|)  »  •„  -  X  f  F ,  (13) 

-e(X)  ■  In  p(X  I  (<k)E*|)  <  «„  -  *  <  F  .  04) 

Based  on  the  log  likelihood  decision  function,  the  feasible  track  set  is 
F  •  j  X|X  Zm  In  p(X  I  { &k )  j|B| )  >  | ■  US) 

Primarily,  the  only  random  component  of  c(  X )  is 

S4'V*'4  .  (16) 


which  for  real  tracks  Is  a  chi-squared  random  variable  with  n  -  dimension 
(1)  degrees  of  freedom.*  Therefore,  error  probabilities  can  easily  be  com¬ 
puted  for  the  hypothesis  test  to  predict  the  accuracy  of  feasible  track  con¬ 
struction.  This  has  a  critical  impact  on  the  ultimate  accuracy  of  the  track¬ 
ing  algorithm,  since  a  real  track  mistakenly  cxdused  from  F  cannot  be 
used  in  the  subsequent  Bayesian  decision  proces. 


U(S)-0,  S<0 

r(S)  »  e‘s/2U(S). 

U(S)»  1,  s  »  0 

(25) 

where  E  {•  ]  indicates  the  expected  value  operator,  0  Sk  is  the  standard 
deviation  of  Sk,  and  f(S)  is  tire  appropriate  density  function  for  a  two- 
dimensional  random  variable  S .  Define  the  cumulative  chi-squared  score 
as 

"i 

1 

1 

VI 

Wi 

■ 

II 

(26) 

For  easy  evaluation  on  a  single  track  an  evaluation  cost  that  would  be 
quite  meaningful  would  be 

r.  0  — L  s' 

S  Nj  * 

(27) 

For  this  cost  function  the  statistical  parameters  are 

E(Cj)  -  2 

(2«> 

E  ((Cj-2)2)  -  4/Nj  "Cj  '  2/Nj* 

(29) 

<C M-TT1-  y^'^UW- 

2  T(Nj) 

(30) 

For  display  purposes  the  use  of  Cj  as  a  value  measure  of  a  single  track 
makes  a  great  deal  of  intuitive  sense. 

Values  of  Cj «  2  for  reasonable  length  tracks  tend  to  imply  that  the 
filter  parameters  are  set  too  loose.  Tlius,  hy  reducing  Ok  and/or  Rk,  one 
could  obtain  lighter  tracks.  Tracks  Tor  which  Cj  »  2  clearly  represent 
bad  data  assignment  which  should  not  be  kept.  More  precisely,  if 


'Due  to  Hie  random  nature  of  lire  measurement  artiest  time  and  of  the  teaser 
•htth  urates  "the  near  measurement,"  the  development  of  Vk  Is  a  Wo  random 
in  nature  but  M  r»  herd  to  comps'*  lit  nwsnm»  Item  one  nra lira! ton  to  the  neat. 
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C,  >  !♦  «2/N|v* 


when  •  ■  3 (for  a 3-tig  ms  ax)  indicates  either  that  incorrect  data  has 
been  assigned  to  the  track  or  that  tins  Tiller  parameter*  (Q,R)  are  too 
Mia II.  While 


Cj  <  2-o2/N*.  Nj  >  3  for  •  •  3 


indicates  that  the  filter  parameters  (Q.R)  are  set  too  large. 

For  a  total  surveillance  region  picture  of  M  measurement  points  and 
L  tracks,  one  has: _ 

j  Nwokrr  of  Psion 

Track  No. _ Scare  in  Track _ 

i  ci  Si 


Track  No. 

Score 

1 

2 

C2 

L 

CL 

The  score  for  the  total  area  should  be  made  up  of  these  scores  and  the 
cost  for  assigning  a  point  to  clutter.  The  number  of  points  in  clutter  is 
Nq,  the  number  of  points  assigned  to  targets  is  Np,  and  the  total  number 
of  measurements  are  M.  These  are  related  by 

L 

M  «  Nc  +  Nj>  Np  *  £  Ni  ■  (33> 

i-l 

A  meaningful  score  could  be  defined  as: 

L 

S*£NiCj*NcSc  (34) 

i-l 

where  Sc  is  a  score  defined  for  a  clutter  point.  If  all  the  measurements 
are  correctly  assigned  to  the  track: 

L 

E(S)  ■  2  £  Nj  +  N(;sc  (35) 

i”l 

»  2(M  -  Nc)  +  NcSc .  (36)- 

If  we  define  the  clutter  score.  S^-,  and  a  total  surveillance  score  as 
C  -  ±  S  Of) 

M 

where,  if  all  the  measurement  points  belong  to  the  track, 

E(C)  -  1  E(S)  *  tt  1 2(M  -  Nc)  ♦  Nc  Scl  (38) 

M  M 


E(C)  -  2 


M-NC  ,  Nc  ; 


so  that  a  plot  of  the  score  for  a  global  surveillance  region  for  s  teal  case 
of  Np  measurements  assigned  to  N*  targets  with  all  points  assigned  cor¬ 
rectly  and  Nc  clutter  points  can  be  geometrically  described  by  Figure  I. 

The  curve  from  N«0,  that  is,  all  points  assigned  to  clutter,  to  N"Np, 
the  correct  number  of  points  assigned  to  clutter,  is  just  a  straight  line 
described  by: 


Ng  varies  from  M  to  M-Np . 

Beyond  the  point  (N“Np),  a  clutter  point  must  be  assigned  to  the  track 
(asMiming  just  one  track  in  clutter).  The  cost  to  assign  a  single  clutter 
point  to  the  track  can  he  calculated  in  lire  following  manner:  A  given 


track  wilt  project  to  a  given  iuhiiI  in  mcusutcttH-nl  space.  I-im  rsatiipk . 
consider  the  two-dimensional  measurement  vector  ( r .  •  I  described  in 
Figure  2- 

If  a  chiller  point  is  added  to  the  track,  the  increase  in  score  on  the 
track  will  be  given  by  Equations  (2 1 )  and  (37).  The  dtttcrence  is  tlia< 
now  the  value  of  the  score  can  be  written  as 

(  L  ) 

(  i»l  1 

uwac  all  track 
points  arc  properly 
Mutual 

Here  Cf  is  the  cost  for  assigning  one  clutter  point  to  the  track. 
is  a  random  variable  and  its  distribution  depends  on  the  distribution 
of  the  clutter  points.  From  Equation  (21 ). 

Cl  "  Uk-lkJV-k'Uk-^)  (42) 

where  z£  are  the  r.p  points  for  the  clutter  point.  The  distribution  of 
this  random  variable  depends  on  the  random  nature  of  the  duller. 

If  the  distance  of  the  duller  point  from  the  predicted  point  a  assumed 
to  be  Gaussian  in  r  and  •  with  aero  mean,  the  distance  from  the  pre¬ 
dicted  point  wilt  be  Rayleigh  distributed  and  the  cost  or  score  C|  ,  the 
weighted  square  of  the  distance  (42),  will  have  an  exponential  density 

1 W._l_  ,.«•=„ M;  “‘'’•“•"‘J  ,41, 

C1  v02  u(c)  “  1.00 

where  o-  is  a  measure  of  the  dispersion  of  the  clutter  points  with 
respect  to  the  mid-  or  trackcr-prcdicted  point.  The  uncertainly,  o, 
takes  into  account  the  unequal  variance  in  r  and  •  and  is  given  by: 

2  (range  of  o  in  surface)  (range  of  r  in  surface) 

®  (t  il 

Vt 

The  score  for  the  closest  clutter  point  to  a  point  predicted  by  the 
tracker  (closest  in  the  sense  of  having  the  smallest  score  defined  by 
(42))  will  be  distributed  as  follows  (x  »  score  of  the  closest  point): 

F„(x)  *  probability  that  the  closest  clutter  point  out  of 
having  a  score  less  than  or  equal  to  X. 

Fx(x)  ■  (4$) 

From  the  theory  of  extremals,  it  is  obvious  tliat  this  distribu¬ 
tion  is  identical  to  the  probability  distribution  of  the  event  defined 
below, 

Fx(x)  »  {  probability  that  all  or  the  Nc  clutter  points  have 
scores  less  than  or  equal  to  x) 


F„(x)  ■  I  -  {probability  that  all  Nc  clutter  points  have 
scores  greater  than  x) 

From  (43),  it  follows  that  the  probability  that  the  score  of  any  one  of 
the  Nc  clutter  points  being  greater  than  x  is  given  by: 

f-x/2oJ  .  (46) 

Assuming  the  clutter  points  are  independent  of  one  another,  the 
probability  that  all  of  the  Nc  clutter  points  have  scores  greater  than 
x  is  given  by 

t-Ncx/2o3 

and  the  probability  that  alt  of  these  clutter  points  have  scores  less  than 
or  equal  to  x  is: 

•Ncx/2o2  .... 
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The  expected  value  of  the  cost  of  aligning  the  closest  clutter  point  to  a 
track  it  then  given  by: 


E(x) 


■/; 


xf(x)  dx  mlSL^ 
Nc 


(49) 


The  expected  coat  for  assigning  the  next  nearest  point  to  a  track  is 

?a2 


(N'C  •) 


(SO) 


and.  finally,  the  expected  cost  for  adding  the  k,h  closest  clutter  point 
out  of  N£  to  a  track  is  given  by: 


Nc*  l-k 


(SI) 


Note  that  in  computing  the  expected  value  of  the  total  cost  function 
defined  in  (41 ).  one  must  also  account  for  the  linear  decrease  in  cost 
caused  by  the  decreasing  weighting  coefficient  (Nc-k)  on  Sc  when  k 
clutter  points  are  incorrectly  assigned  to  tracks.  This  amounts  to  a 
decrease  of  Sc/M  for  every  clutter  point  assigned  to  a  track.  Therefore, 
the  net  increase  in  the  expected  cost  by  assigning  the  k  closest  clutter 
points  out  of  Nc  to  tracks  is  given  by : 


In  clleel,  this  Is  Hie  sum  ol  an  hivicesing  liypeilMtoc  IuikIkmi  amt 
a  decreasing  linear  function. 

Stowe,  on  the  average,  »Sc.  the  shape  of  the  right-hand  aide  of 
the  curve  in  Figure  I  is  concave  and  mo  notonic  ally  increasing. 

Note  the  two  regions  of  this  figure.  In  one  part  fewer  points  arc 
aesified  to  tracks  than  actually  arc  available,  i4.,  loo  many  detections 
were  missed.  In  the  other  half  of  the  figure,  false  ataims  or  chitter 
points  were  assigned  to  tracks.  In  both  pans  of  the  curve  the  assump¬ 
tion  is  made  that  the  data  association  is  done  in  an  optimal  manner  on 
the  average.  Thus,  if  a  clutter  point  is  added  to  a  track,  it  is  the  clutter 
point  that  lies  “closest''  to  a  track  of  those  unaasigned.  The  curve  also 
assumes  that  all  measurements  correctly  assigned  as  track  points  are 
assigned  to  the  correct  track.  Wrong  assignments  of  data  would,  of 
course,  make  for  even  worse  scores  “on  the  average.” 

4.  CONCLUSIONS 

In  the  last  section,  it  was  shown  that  the  refined  scoring  algorithm 
possesses  appealing  properties  from  a  geometric  viewpoint.  There 
exists  a  unique  number  of  points  in  a  track  that  results  in  lower  average 
cost  than  any  other  number.  Also,  the  sensitivity  of  the  score  to  varia¬ 
tions  in  assigned  number  of  track  points  can  be  controlled  by  the  duller 
score,  S£.  This  is  readily  apparent  from  Figure  I. 

This  scoring  algorithm  is  therefore  a  very  useful  approach  in  extract¬ 
ing  clutter  points  out  of  a  given  target  track. 

The  algorithm  is  currently  being  applied  to  an  ocean  surveillance 
problem  and  the  results  of  this  are  very  encouraging.  For  a  given  data 
set  having  a  false  alarm  rate  of  1CT4.  i.e.,  one  clutter  point  in  every  I04 
measurements,  using  the  refined  scoring  algorithm  defined  in  Section  3. 
we  have  found  tliat  we  can  effectively  decrease  this  false  alarm  rate  to 
KT7.  This  represents  a  three-order-of-magnitude  decrease  and,  hence, 
the  detection  capabilities  of  the  tracking  algorithm  have  been  signifi¬ 
cantly  enhanced. 


> 


N  •  NUMgfn  OP  POINTS  ASSIGNID  TO  TtoACKS 


figure  I.  Geometries  description  of  refined  tearing  algorithm. 
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This  paper  addresses  the  laaue  of  target  tracking 
when  confronted  with  a  aet  of  aound  apeed  parameter! 
that  are  partially  or  completely  unknown.  It  explore* 
the  caae  where  theae  paraeeter*  are  augmented  to  tha 
target  atate  In  an  extended  Kalman  filter.  The  filter 
proceaae*  measurements  of  aound  tlne-of-arrival  differ¬ 
ence  and  Doppler  difference  f roe  a  aet  of  apatlally 
dlaplaced  eenaora. 

For  acenarlo*  Involving  up  to  three  aenaora  it  haa 
been  found  thet  blaaed  target  poaltlon  estleatea  and 
narglnal  lyataa  obacrvablllty  occur*.  Thla  1*  readily 
verified  by  propagating  tha  elgenvaluea  of  the  infor¬ 
mation  matrix  In  tine.  Using  thla  aa  an  analyil*  tool, 
a  number  of  geoeatrlcal  aenaor  configuration*  are 
analyxed. 

In  general.  It  la  found  that  with  three  aeaiora,  the 
system  la,  at  beat,  marginally  obaervable  for  any  geom- 
etry.  However,  when  ualng  four  and  nor*  aenaora,  ays- 
ten  obacrvablllty  and  estimation  performance  are  narkedly 
Improved  when  two  of  the  aound  speed  filter  paraeeter* 
arc  specified  to  within  a  close  tolerance  of  their 
actual  values.  When  attempting  to  estimate  all  of  the 
sound  speeds  (or  for  that  matter  (n-1)  sound  speeds, 
n  •  number  of  sensors) ,  It  Is  again  noted,  at  in  the 
three-sensor  case,  that  system  observability  and  estima¬ 
tion  performance  become  degraded. 


1.  Introduction 


lo  a  target  tracking  application,  one  reason  for  poor 
estimation  performance  can  be  a  lack  of  knowledge  con¬ 
cerning  the  parameters  of  the  mathematical  model  that 
relates  the  target  state  to  the  measurements.  Since  a 
mathematical  model  Is  a  necessary  ingredient  to  any  tar¬ 
get  cracker  or  state  estimator,  use  of  incorrect  parame¬ 
ters  could  lead  to  estimates  that  diverge  from  "truth" 
over  a  period  of  time. 

It  is  this  problem  that  Is  dealt  vlth  In  this  psper. 
More  specifically.  It  Involves  a  target  tracking  prob¬ 
lem  where  measurements  of  time  difference  and  Doppler 
difference  are  collected  from  pairs  of  spatially  dis¬ 
placed  sensors.  The  central  Issue  Is  that  the  sound 
speeds  from  the  target  to  each  of  theae  sensors  are 
partially  or  completely  unknown.  These  speed  psrameters 
appear  in  the  mathematical  model  relating  the  target 
state  to  the  measurements.  To  avoid  biased  target 
tracks,  some  mechanism  should  be  found  to  accosmodat* 
these  parameter  uncertainties. 

The  material  In  this  paper  Is  basically  an  extension 
of  an  earlier  work  ( 1 )  and  is  more  conclusive  In  terms 
of  the  results  that  were  obtained  for  a  number  of  dif¬ 
ferent  geometric  scenarios  Involving  ssnsor  placements 
and  target  location. 

Ue  will  start  our  discussion  In  Section  2  by  describ¬ 
ing  the  mathematical  model  for  the  given  process  and 
proceed  to  define  an  estimator  that  can  accommodate 
both  unknown  sound  speeds  and  the  target  state  vector. 
All  of  the  simulation  results  will  be  presented  In 
Section  3.  In  addition,  we  will  also  show  how  we  can 
assess  system  observability  via  the  information  matrix. 
This  will  play  a  useful  role  in  exploring  system  observ¬ 
ability  for  a  number  of  different  target/sensor  geome¬ 
tric  scenario*. 


From  the  simulation  results  It  will  ba  scan  that  a 
minimum  number  of  aenaora  are  needed  and,  la  addition, 
two  of  the  aound  speeds  must  be  known  correctly  before 
accurate  estimation  of  the  target  state  can  be  achieved 
Finally,  In  Section  A,  a  a  usury  of  the  result*  and 
conclusions  will  be  presented. 


2.  System  Definition 


One  area  Where  the  uncertain  model  parameter  problem 
could  arise  is  depleted  In  Figure  1.  It  could  equally 
well  apply  to  tracking  of  vehicles  on  the  Earth  via 
geophonea  or  any  place  that  the  signal  does  not  travel 
vlth  an  Infinite  effective  or  known  velocity.  Ue  have 
a  set  of  1  spatially  displaced  sensors  (distance  to 
target  *  R^) .  If  the  target  generates  or  reflects  sound 
at  time  Instant,  t,  because  of  the  sound  travel  time  to 
each  sensor  (sound  speed  -  e*),  It  will  be  seneed  at 
each  sensor  at  times  tj,t2,...,ti.  In  addition.  If  the 
target  saves  at  a  velocity,  v,  the  Doppler  sensed  at 
each  sensor  will  be  different.  To  estimate  the  target 
state,  l.e.,  position,  speed,  and  course,  measurements 
of  sound  time-of -arrival  difference  and  Doppler  differ¬ 
ence  for  a  sensor  pair  1-J  (1,3*1, 2, 3, . . .  ,1*3)  can  be 
processed  through  a  Kalman  filter.  Using  spherical 
geometry,  these  measurements  can  be  related  to  the  tar¬ 
get  state  by  the  following  equations  (!]: 


'  ■  t.-t.  - 

Rt  ^ 

(1) 

13  1  3 

cl'c3 
/  R.  R.  \ 

'  •  f  -f  - 

(--U-l)  f 

(2) 

13  1  3 

\  c*  CW 

where 

R i  -  cos  ^slnXjSinX^^  +  cosXjCosXjCosfXj-B^) ]  (3) 

XjlsinXjCO*XJcos(x2-ei)  -  cosxjSlnX^] 

R1  “  =7  «  sin  Rt 


x4eosx1cosX1sin(Xj-91) 
sin  Rt 


(A) 


In  the  above  equations,  xj., 1-1 ,2, 3, A  represent*  the 
target  latitude,  longitude,  latitude  rate,  and  longitude 
rate,  respectively.  The  latitude  and  longitude  of 
hydrophone  1  Is  Xj  and  8j.  An  Implicit  assumption  In 
(A)  Is  that  the  sensors  are  stationary. 

Throughout  the  paper  ve  will  take  the  target  atate 
vector  to  be  1(7  •  Ixj .x2.x3.x4] ,  The  rceson  for  doing 
this  is  that  target  motion  can  be  described  by  a  linear 
set  of  equations.  In  discrete-tlxM  form,  the  equation* 
for  the  target  dynamics  are  given  by: 
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where  At  *  eampllng  Interval;  l.a.,  the  tine  between 
ncasurenents  of  Tij .fjj ;  w(k)  is  e  zero-naan,  white 
noise  veetor  sequence  that  perturbs  the  target  fron 
otherwise  constant  courae/speed  notion;  and  where 


E(5(k)  ^(j)}  -  Qk  6^  . 


1c  is  a  staple  natter  to  display  speed  end  course  at  any 
tine  by  using  the  equations  below: 


COUISE  -  ten 


Since  the  naaaureaent  nodal  is  nonlinear  (eq.  1-4), 
one  can  laplenent  an  extended  Kalnen  filter  to  track  or 
obtain  estlnates  of  the  target  state  vector,  x(k).  This 
is  easily  done  by  linearizing  the  aeeaurenent  equations 
about  the  aost  current  state  estimate,  x(k),  to  obtain  a 
linear  neasurenent  equation,  H(x(k)).  The  equations  for 
the  filter  are  standard  (2)  and  are  summarized  below: 


2(k+I/k) 


♦(At)  x(k/k) 


P(k+l/k)  -  d(At)  P(k/k)  9*(At)  +  Qk 
x(k+l/k+l)  -  x(k+l/k)  +  ^jlztk+l) 
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P(k+l/k+l)  -  11  -  k,H(x(k+l/k)))  P(k+l/k) 
-  P(k+l/k)  H(x(k+l/k))  [H(2(k+l/k)> 
P(k+l/k)  HT(x(k+l/k))  + 


z  ;k)  -  [t13 Ck),  f13 (k)) 
hT(S(k+l/k))  -  It^tk).  fi3(k)]  - 
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R,  -  covariance  matrix  of  additive  white 
noise, that  contaminates  z(k) 


H(£(k+l/k)) 


3h(T(k)l 

3x(k) 


T(k)  -  5(k+l/k) 


3.  Simulation  Reaulte 

To  cxaalne  the  effects  of  estlaatlng  unknown  sound 
speeds,  we  selected  a  number  of  different  cases  involv¬ 
ing  different  target  notion  scenerlos  and  three-sensor 
configurations  as  shown  in  Figure  2. 

The  locations  of  the  sensors  are  defined  in  Tabic  1 
and  the  target  motion  scenarios  (cases  1-24)  are  sum- 
narlzed  in  Table  2. 

Table  1.  Location  of  sensors. 


Array 

Latitude,  X 

Longitude,  e 

1 

3  deg 

0  deg 

2 

-2.5  deg 

4.33  deg 

3 

-2.5  deg 

-4.33  deg 

Table  2. 

Target  notion  scenarloa. 

Case 

Humber 


From  (1),  (2),  (13),  and  (14)  it  is  easy  to  see  how 
sound  speed,  eg,  enters  into  the  measurement  model.  In 
(1],  it  wee  shown  that  when  incorrect  values  of  Cg  were 
used  in  the  filter  model  (assuming  one  did  not  know  the 
true  c^),  the  resulting  state  estimates  were  found  to  be 
biased  off  from  the  true  target  state.  In  some  cases 
these  biases  were  significant,  and  consequently  the 
deviation  fron  truth  was  as  significant. 

To  compensate  for  this  problem,  the  sound  speeds  vere 
treated  as  additional  state  variables  and  augmented  to 
the  target  state.  Since  the  sound  speeds  were  constant 
over  the  estimation  Interval,  the  state  dynamics  were 
simply  defined  by  tg  •  0.  The  extended  Kalman  filter 
was  then  implemented  for  this  augmented  state  veetor  to 
generate  estimates  of  both  the  targec  state  and  the 
unknown  sound  speeds. 

In  the  next  section,  we  will  summarize  some  of  the 
earlier  results  that  were  obtained  and  then  present  more 
exhaustive  results  that  indicate  a  definite  trend 
occurring. 


Latitude 
■  3  deg 


Longitude 
.3  deg 


Speed 
10  knota 


2.5  deg 


-1.3  deg 


T 

-1.3  deg 


T 

-2.0  deg 


Couree 
0  deg 
90  deg 
ISO  deg 
270  deg 
0  deg 
90  deg 
ISO  deg 
270  deg 
0  deg 
90  deg 
180  deg 
270  deg 
0  deg 
90  deg 
180  deg 
270  deg 
0  deg 
90  deg 
180  deg 
270  deg 
0  deg 
90  deg 
180  deg 
270  deg 


He  made  the  following  assumptions: 

(a)  The  covariance  matrix  of  the  discrete-tine 
process  dynamics  was  defined  by: 


a2  M. 

q44  2 


T 

-4  deg 


Q(k)  - 


2  At* 
q33  2 


q44  6t 


V  * 


where  qjj  •  q44  •  (.0091827) 2  knece/eee  Is  the  power 
spectral  density  of  the  rendoa  nolee  perturb ing  the 
velocity  stete  equations  of  the  continuous  eystea.  This 
randowess  In  the  target  velocity  for  the  continuous 
systea  translates  both  Into  s  position  and  velocity 
uncertainty  In  the  equivalent  dlacrete-tlae  model.  The 
values  of  qjj  and  q44  roughly  correspond  to  a  seandsrd 
deviation  of  .34  nautical  alias  In  position  and  .6  knot 
In  velocity  over  a  time  interval  of  one  hour  In  the 
dlscrete-tlae  aodel. 

(b)  The  aeasureaent  matrix  covariance  matrix  was 
defined  by: 


(c)  8t-300  sec  was  the  nominal  time  Interval  between 
aeasureaents . 

(d)  The  filter  processed  measurement!  froa  the  sensor 

pairs  In  a  sequential  Banner  starting  with  sensor  pair 
1-2.  1-3,  2-3.  1-2,  1-3 . etc. 

(e)  The  sound  speeds  from  the  target  to  each  of  the 
tensors  were  chosen  as  [3]: 

Cj  -  4857  ft/aec 
c2  ■  4850  ft/sec 
Cj  -  4870  ft/sec 

Ve  started  out  by  assuming  that:  first,  only  one  of  the 
three  sound  speeds  was  unknown  and  consequently  was  esti¬ 
mated  along  with  the  targec  state;  second,  two  sound 
speeds  were  unknown  and  were  estimated  along  with  the 
target  state;  and  third,  all  three  sound  speeds  were 
unknown  and  estimated  along  with  the  target  state.  In 
all  of  these  cases.  It  was  found  that  biased  estimates 
were  generated  by  the  tracker.  A  typical  example  of  this 
is  shown  In  Figures  3,  4,  and  5  where  an  attempt  was  aade 
to  estimate  the  unknown  sound  speeds  snd  the  target  state. 
The  solid  curves  represent  the  truth  model  whereas  the 
dashed  curves  represent  the  state  estimates.  Note  the 
significant  biases  In  latitude  and  longitude  in  two  of 
the  sound  speeds. 

Because  of  these  biases,  It  was  decided  to  examine 
the  observability  of  the  system  for  all  of  the  cases 
defined  in  Table  2.  This  is  easily  done  with  the  aid 
of  the  lnformatloh  matrix  1 2).  For  the  case  involving 
no  process  noise  and  state  vector  a  priori  Information, 
the  information  matrix  Is  Identical  to  the  Inverse  of 
Che  Kalman  filter  covariance  matrix,  P-1(k/k).  This 
matrix  must  be  positive  definite  for  stochastic  observ¬ 
ability  and,  provided  the  above  conditions  apply,  is 
given  in  recursive  fora  by: 

F_1(k/k)  -  ♦Vit)  P*"1  (k-l/k-1)  (-at) 

+  HT(x(k-l))  R^HCxlk-l)); 
p"'(0/0)  -  0  (15) 

where  *(&t)  is  the  state  transition  natrlx  defined  In 
(5),  HOf(k-l))  Is  the  measurement  matrix  linearized 
about  the  stats  vector  x(k-l). 

To  assess  the  property  of  stochastic  observability, 
the  normalized  eigenvalues  of  this  matrix  ware  computed 
(normalized  to  one)  and  plotted  as  a  function  of  time. 
Figure  6  shows  the  results  that  were  obtained  for  case 
11  in  Table  2.  One  of  the  position  eigenvalues  becomes 
ill-conditioned  and  exhibits  s  smaller  maximum  magni¬ 
tude  than  the  other  position  eigenvalue  by  a  couple 
orders  of  magnitude.  This  analysis  was  repeated  for  ell 
of  the  other  23  cases  and  the  same  general  result  was 
obtained,  l.e..  Ill-conditioned  behavior  of  one  of  the 
eigenvalues.  Because  of  this  and  the  fact  that  the 
state  estimates  were  blesed,  we  concluded  that  the 


system  was  marginally  observable  for  a  three-sensor 
configuration  and  unknown  sound  speeds. 

As  s  means  of  enhancing  system  observability,  it 
was  decided  to  Introduce  more  chan  three  sensors  for 
tsrget  tracking. 

Ua  first  started  with  four  sensors  using  different 
sensor/target  motion  geometries.  Four  cases  were  con¬ 
sidered  and  the  geometries  are  aiunarized  In  Figures 
7  to  10. 

Using  the  same  philosophy  as  In  the  three sensor  case 
earlier,  we  started  out  by  estimating  one,  two,  three, 
and  then  four  sound  speeds.  For  all  of  these  cases. 

It  was  found  that  we  could  estimate  the  target  etate 
and  up  to  two  sound  speeds  without  obtaining  biased 
estimates,  but  as  soon  as  we  attempted  to  estimate 
three  or  four  sound  speeds,  biases  in  the  estimates 
again  were  noted.  Marginal  system  observability  again 
was  suspect.  To  aubstantlate  this  we  looked  at  the 
eigenvalues  of  the  information  natrlx  as  a  function  of 
time.  The  functional  variations  of  the  eigenvalues  were 
found  to  be  relatively  smooth  and  monotonlcally  Increas¬ 
ing  for  estimation  of  one  or  two  sound  speeds.  An  exaar' 
pie  of  this  is  presented  In  Figure  11.  It  Involved  the 
target/hydrophone  geometry  defined  by  Figure  10  where 
we  estimated  the  target  state  and  two  of  the  sound 
speeds.  However,  as  we  began  to  estimate  three  end 
more  sound  speeds,  the  function  variation  of  several  of 
the  information  matrix  eigenvalues  becomes  progressively 
more  ill-conditioned  and  lower  in  absolute  magnitude— 
an  Indication  that  the  property  of  system  observability 
has  been  weakened. 

To  complete  our  analysis,  we  then  explored  the  use 
of  five  sensors.  Two  geometries  were  selected  and  are 
shown  in  Figures  12  and  13. 

Using  the  same  approach  as  before,  we  began  by  esti¬ 
mating,  first,  one  sound  speed,  then  two  sound  speeds, 
and  so  on.  Interestingly  enough,  it  was  found  that  one 
could  now  estimate  up  to  three  sound  speeds  before 
biased  estimates  again  occurred. 

For  all  of  the  above  cases  Involving  four  and  five 
receiving  sensors,  the  general  observation  was  that  one 
ceuld  estimate  the  target  state  and  up  to  two  sound 
speeds  for  the  four-sensor  configuration,  and  the  target 
state  and  up  to  three  sound  speeds  for  the  five-sensor 
configurations. 

Of  course,  these  conclusions  are  based  upon  a  finite 
set  of  examples,  and  to  substantiate  the  above  claim 
more  rigorously,  one  would  have  to  Implement  a  more 
exhaustive  sec  of  examples. 


4,  Conclusions 

In  summary.  It  was  first  noted  that  target  tracking 
via  extended  Kalman  filtering  tends  to  produce  biased 
estimates  when  the  sound  speeds  were  uncertain  and 
Incorrectly  specified  in  the  filter.  Attempts  to  Addi¬ 
tionally  estimate  the  sound  speeds  were  shown  to  be  of 
no  avail  In  eliminating  these  biases— even  when  Apply¬ 
ing  traditional  filter  parameter  variations  that  In  past 
applications  tended  to  make  the  filter  more  robust  to 
parameter  uncertainties. 

For  this  reason  the  observability  of  the  system  was 
explored  in  greater  detail.  With  the  aid  of  the  Informa¬ 
tion  matrix.  It  was  found  that  the  system  was  marginally 
observable  over  the  geographical  region  defined  by  the 
three  receiving  sensors. 

Because  of  this,  we  therefore  took  a  look  at  using 
time-difference  and  Doppler  difference  measurements 
froa  more  than  three  sensors.  In  particular  we  looked 
sc  configurations  Involving  four  and  live  receiving 
sensors. 

The  resulcs  froa  a  finite  set  of  examples  have  shown 
that  target  tracking  performance  Is  Improved,  l.e.,  very 
small  or  nonexistent  biases,  but  estimation  of  all  sound 
spaads  Is  not  possible.  Cenerally  speaking.  It  seems 
that  If  we  were  given  n>3  receiving  arrays.  It  would 


b«  poeelble  to  ncUili  ch«  target  atata  and,  at  Boat, 
n-2  of  tha  aound  apaada  to  aaeh  of  ehaaa  aanaora .  Tha 
remaining  two  aound  apeada  have  to  ba  apaelftad  a  priori 
for  the  flltor. 
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Figure  2.  Target  notion  acenarloa  for  aound  epeed 
eetlaatlon. 
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Figure  1.  Target  atata  and  aound  apeed  eetlaatlon 
(latitude  and  longitude). 
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Figure  9.  Geoaecry  of  four  censor*-- Cat*  3. 
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Figure  10.  Ceonetry  of  four  sensors— Casa  4 
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Figure  11.  information  ncrli  eigenvalues— target 
state  and  two  sound  speeds. 
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Figure  12.  Ceosetry  of  five  tensors— Case  1. 


Figure  I).  Ceooetry  of  five  eensort— Case  2 
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Summary 

The  concept  of  coherence  ea  a  maaaura  of  the 
linear  relationahip  between  two  or  more  time  ranee 
ia  diacureed.  The  definition  of  coherence  in  term*  of 
the  complen  eroee  power  apectral  denrity  matrix  ie 
given  and  the  relationahip  between  pairwlaa  and  multi* 
pie  coherence  and  channel  eignal  to  noiae  ratioa  ia 
diacuaaed.  A  aamplc  atutiatic  for  coherence  te  given 
while  the  detaila  daacribing  tha  performance  which 
can  be  obtained  from  thie  atatiatic  are  contained  in 
a  aeparate  report  by  theae  authora1. 


1.  Introduction 

A  problem  of  lntereat  in  many  different  dice* 
tplinea  ia  that  of  determining  if  there  ie  a  met  curable 
relationahip  (phyaical  cauarlity)  between  two  er  more 
lime  aeriea,  in  addition,  one  would  often  like  to  ob¬ 
tain  a  quantitative  meaningful  meaaure  of  the  degree 
of  that  relationahip.  Thia  paper  doacribea  one  poaai- 
hie  meaaure  of  auch  a  relationahip,  the  coherence. 

The  moat  common  meaaure  of  auch  a  relation¬ 
ahip  ia  the  pairwiae  or  multiple  correlation  coefficient. 
The  nature  of  the  correlation  coefficient  ia  well  docu¬ 
mented  and  will  not  be  diecueecd  here  other  than  to 
note  that  tt  ia  not  a  function  of  frequency  and  may  be 
affected  by  linear  tranaformatione  of  either  of  the 
time  aeriea. 

Tne  coherence  function  (magnitude-equared 
multiple  or  pairwiae  coherence  function)  ia  defined  aa 
a  frequency-dependent  quantity  that  rangca  between 
aero  and  one**1  .  Thia  coherence  function  ia  aero  if 
the  two  or  more  ergodic  time  aeriea  are  independent 
tuncorrelated,  if  Cauaaian)  and  equal  to  one  at  any 
frequency  where  there  ia  a  linear  tranaformation 
between  the  one  or  more  input  time  aeriea  and  the 
output  or  reference  time  aeriea. 

The  aituatton  of  intereat  ia  ehown  in  Figure  1 
where  V,  indicatee  the  noiae  contaminating  the  algnel 
u(ii  in  the  Ith  channel.  In  general,  each  tranamiaaion 
channel  ia  compoaed  of  linear  and  nonlinear  parta 
(Figure  2).  The  eum  of  the  output  of  the  nonlinear 
ayatem  (uaually  a  email  part  of  tne  total  tranamiaaion), 
tne  meaeurement  noiae,  and  the  background  noiae  ia 
.rouped  into  the  effective  noiae  term  Oi(t)  (Figure  3). 

We  are  tntcrcated  in  detecting  the  preaence  of 
a  common  aignal  u(t)  tn  two  or  more  chatuiela.  Tha 
mpai-output  relationahip  indicated  in  Figure  1,  3  can 
oe  writier,  aa 


a 

«(t)  •  / 

-•fi i 


«dr  ♦  *(»). 


Note  that  becauoe  of  phyaical  cauaality  require- 
monte  g(r)  ia  aoro  for  all  r  leaa  than  aero,  lb  fact,  it 
will  be  aero  for  all  r  laaa  than  aoma  poeitive  time 
which  ta  tha  time  it  takoa  a  aignal  to  travel  from  the 
aource  to  a  eenaor, 

Tha  true  value  of  the  coherence  between  time 
aeriea  ia  generally  an  unknown  quantity.  In  fact,  any 
meaaure  of  the  relationahip  between  two  or  more  time 
aeriea  generally  muet  be  baaed  on  time  tracaa  of  thoaa 
aeriea.  The  functional  relationahip  between  the  time 
aeriea  and  the  meaaure  or  eattmate  of  coherence  ia  a 
aamplc  atatiatic  for  coherence.  The  aoeumption  that 
any  one  infinite  length  aampla  of  each  aeriea  will  be 
enough  to  allow  ue  to  eetimate  the  coherence  exactly 
ia  made  implicitly,  and  all  time  aeriea  are  aaaumod 
to  be  atationary  and  ergodte.  Unfortunately,  in 
practice,  one  ia  given  only  a  finite  amount  of  data 
from  each  of  the  time  aeriea,  In  thia  caae  the  aamplc 
atatiatic  ie  a  random  variable  dietributed  about  the 
"true  magnitude-equared  coherence".  The  denaity 
function  for  the  aamplc  atatiatic  in  thia  report  ia  given 
for  a  number  of  apecific  valuea  of  degreea  of  freedom 
(N)  and  number  of  time  aeriea  (M)  in  reference  1. 

Baaed  on  theae  denaity  functiona  receiver 
operation  cnarectenatic  (ROC)  curvea  have  been  de¬ 
veloped.  Theae  curvea  define  the  probability  of 
detection  veraua  probability  of  falac  alarm  for  a 
aignal  of  a  given  true  coherence.  Curvea  of  probabil¬ 
ity  of  detection  veraua  true  coherence  for  fixed  levela 
of  probability  of  falac  alarm  have  aleo  been  developed 
for  many  degreea  o!  freedom  and  number  of  eenaora 
up  to  ID.  Theae  are  ail  reported  in  detail  by  theae 
eulhorc  in  reference  1.  The  fundamental  queation  of 
the  relationahip  between  input  eignal  to  noiae  levela 
and  true  coherence  for  the  two  channel  and  multi¬ 
channel  ceaea  of  Figure  1  ia  diacuaaed  in  thia  paper. 

For  a  general  diecuaaion  of  the  concapt  of 
coherence,  the  reader  ia  directed  to  reference  3. 

For  a  detailed  derivation  of  tha  diatribution  of  the 
patrwtee  and  multiple  coherence  atatiatic,  the  reader 
ia  referred  to  reference  7,  The  daneitiee  and  derived 
performance  curvea  in  reference  1  are  particularly 
difficult  te  obtain  for  low  coherence  and  high  valuea  of 
N  (number  of  aamplea  of  the  time  aeriea).  and  baaed 
on  the  authora'  knowledge  are  not  available  elaewhere 
in  the  literature. 
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Croaa-Power  Spectral  Density 
Matrix  and  MwIlM*  Cnhnrnnfn 


Tk&a  la  the  mutual  or  pairwiat  coherence  between 
channela  one  ami  two. 


Multiple  coherence  c«n  be  moat  •tally  defined 
m  term*  of  the  croaa  apectral  denaity  matrix 

where  the  element#  of  thia  matrix  are  defined 

Vjy 


<i.») 


Th«  ctoiipowtr  apectral  dtulty  matrix  ia  of 
court*  equivalent  to  tha  croiacorrclatiOB  matrix. 

Eithtr  of  thaaa  togathar  with  tha  mean*  of  tha  M 
jointly  gauatian  a  tat  ionary  procaaaaa,  x  (t),  x,(t), . . , 
x  (t|,  completely  apnei/i .«*  tht  joint  diatribution 
function  of  thaaa  procaaaaa. 

Civan  M  finlt*  longth  tlma  tracoa,  thora  ar* 
wall  known  technique*  /or  obtaining  "aampl*  oatimatoa" 
of  tha  croaa*  and  autopowar  apoctral  alamanta.  That* 
catimatta  ar*  ueed  to  obtaui  aampl*  aatimatta  far  tha 
multipl*  cohtranca  bttwatn  tht  varieua  tim*  aaria*. 

Tha  aampl*  aatimat*  for  tho  croea-power 
apectral  danaity  matrix  ia  a  function  of  the  baaic  data 
x  (t).x  (t). .  ..x  ft)  ovar  aoma  finit*  tim*  record.  It 
data  not  raquiraanowlodg*  of  any  of  the  charactar- 
iatica  of  th*  tranamiaaioo  channel*  or  of  tha  aignal-to- 
noia*  ratio*  of  th*  received  aignala.  Th*  meaning  of 
multiple  coherence  will  b*  diacuaaad  in  term*  of  thoa* 
quantitiea  in  latar  taction*  in  an  attempt  to  illuminate 
the  aubjact.  Haro,  however,  w*  will  define  multiple 
coherence  aimply  in  term*  of  th*  croeepower  apectral 
rtenaity  matrix  and  ita  alamanta. 


In  th*  three-channel  c*e*  it  ia  aaaily  aeea 
(writing  Stj(  w)  at  S^l  that 

I  I  *  S»lS2ll  *S 22 IS1 1 1  •i>e15llSl)Sll1 

*  c  e  c  lc  I* 

SllS2*SJj‘Sll|  *gj| 


Uaing  tha  fact  that  pairwiat  coherence  1*  defined 
by  equation  2.  S  wa  can  writ* 

.2  hi .  2^+  hi .  jf"^*^S12S23^Sl  ^®1 1S22*JS* _ 

|*i,2.s|  - - : — s - «*-7» 


MyL2l  ♦frl.il  *  2  R*(yi2  y23y31^ 

1*3 


Ndttthat  if  the  croaapower  apectral  density  of 
channel*  two  and  three  it  aero  (Sgj*0),  the  coherence 
between  theee  two  channel#  ia  aero  and  the  multiple 
coherence  of  channel  one.  given  two  and  three  ia. 

iyi!2.Jr*iyi.2i2+iyi>ji2  ,2-,» 


The  multiple  coherence  between  x.(t)  and 


x  (t),x2(t)... 
{reference  S) 


. .  *M(l)  iaJd*fin*d  by. 


!>j!1.1...J.I,J*1....nK  1-1/  (s..(ui)S,'i(w)].  (2.2) 


«!<.ert  y'lwl  ia  the  jth  diagonal  element  of  the  inveree 
The  multiple  coherence  of  the  jth 


.a:  tr.e  S  iw). 


eenaor  with  reapect  to  the  other  aenaore  repreaenta 
the  proportion  of  the  variance  (power)  of  eenaor  j  that 
can  be  explained  by  a  linear  combination  of  thg  remain* 
int  aenaora  in  a  minimum  mean  aquart  sense  . 

By  reducing  thia  definition  to  tha  aimple  two- 
channel  ca«e  we  can  write 


S2*  r|9,W 


‘ll'-'-SllW 
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3.  .  Pairwiat  Coherence  and 
Signal. to-Noi*t  Ratio* 

If  x  ^  (t )  and  x2(t)  art  generated  ta  indicated 
in  Figure  3.  equation  1.  1  can  be  written  in  th* 
frequency  domain  by  Fouriar  tranaform  at 

X,  ( w)«G  j  ( <4U(  o>)*V  j  ( Vi  (3.1) 

X2(u).C2(w)U(u)*V2(u)).  (3.2) 

Again  aeaumlng  th*  noia*  term*  Vj  and  V  ar*  inde¬ 
pendent,  th*  croee-  and  eutopower  epectral  deneitiee 
can  be  written  aa 


Sn(iu)a  Gj(ei)  2Su(o»+Svl(iu) 

(3.  3) 

S22(«)+  C2(a»)  2Su(o»)aSv2(*i) 

(3.4) 

|S!2(«)|2.  CjlW)2  C2(W,  2Su(0»2. 

(3.3) 

to  that 


S1 ! I u )»S22(«)’[Sj j (u»)S22(w)  -|S12(w)|1].  (2.4) 


Thia  gtvaa 
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Noun*  that  lc,«tf)|2  5U(W)  and  Ig2(oi)!2  S^eg  or* 
"*»  »I|WI  etenal  power  iptn n  at  iha  receiver  la 

*  1  -*>•  auj  Iwa,  %•  wilt*  llt« 

equared  coherence  a* 


!c,(a»!2  !c2<w!*s  (u»l 
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Da fining  tha  eignal-to-noiae  powar  In  tha  Jth  channel 
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wa  «aa  writa 
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Certain  |tBaril  ittttmnta  concerning  thit 
peirwtee  magnitude -equxrcd  coharaaca,  or  juat  "co¬ 
herence.  "  caa  now  ba  made. 


Tha  coharaaca  ia  bounded  between  taro  and 

one: 

o  *  l>ji2!2  *  » ■  <».») 

1/  tb«  noiaa -to-elgnal  powar  goaa  to  Infinity  In 
aithar  channel,  the  coherence  will  |0  to  aero,  Thia 
will  happen  if  the  aignal  powar  In  that  channel  tadaa 
to  aero.  Alao,  the  naiae-to-eignal  powar  ratio  in 
both  channel,  moat  *o  to  aero  for  the  coharaaca  to  fo 
to  one. 


An  interacting  and  informative  interpretation 
of  the  coherence  between  theae  two  channal  a  can  ba 
made  in  the  following  manner.  Aaaume  one  of  die 
channel,  ia  noiac  free  •  0).  Thia  channal 

then  become,  tha  input  aignal.  The  coherence  be. 
Meet  cr-annela  one  and  two  ia  now  given  by 


'1.2 


iQjfo)!2 

JCj(W)!25u<W)+5vI(Wl 


(S/N) 


1+1S/N),  Sj  +  Nj 


(3.10) 


Thia  magnitude  .aqua  red  coherence  la  tha  fraction  of 
the  power  of  the  output  x,(t)  which  cornea  from  the 
eignel  input  paaeed  through  a  linear  ayatem. 


Since  the  traaemitted  aignal  ia  generally  not 
available,  it  ia  uaeful  to  look  at  thia  phyaical  inter¬ 
pretation  of  coherence  from  a  different  point  of  view. 
Let  ua  eimply  take  aignal  Xjftt  aa  our  baaic  aignal  and 
calculate  the  coherence  between  x,(t)  and  our  given 
a-.gna!  *,(t).  From  equation  2.7,  our  definition,  thia 
will  be  the  aame  reeult  aa  if  we  took  aignal  x,(t)  aa 
our  "given''  aignal.  Thu  a  U(W)  in  equation  3.2  ia  re¬ 
placed  by  XjlU): 

X,  (ell  •  Hj2(0))X2(0))  4  V  ,M)  .  (3.11) 

H,,loi)  ia  tha  effective  linear  tranafar  function  be¬ 
tween  output  Xj(t)  and  output  x,(t).  V^,j(t)  »•  the 


affactive  noiaa  on  tha  traaemieelon  channel.  It  muat 
he  noted  that  ll.JU)  ie  no  longer  neceeearlty  a  caueal 
ayalewi.  Nnw  Ihe  rnherenre  Hetwren  rhennela  one 
and  two  can  be  written  aa  to  aquation  3. 1  Di 
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Thia  two-channel  taagaltuda-aquarad  coherence  ia 
tha  ratio  of  Mm  powar  at  output  x.  (t),  which  ia  canoed 
by  the  "input  x.(t). "  traaemitted  over  the  effective 
linear  traaemiaeioa  channel  [H,-(ei})  to  die  total 
power  in  output  CenalderaDen  of  thia  effective 
linear  traaemiaeioa  channal  allow,  thia  phyaical  in¬ 
terpretation  of  the  coherence  to  be  eaaily  carried 
over  to  multiple  coherence. 


4.  Multiple  Coherence  and 
Simal-to-Moiaa  Ratio, 

In  the  caaa  of  M  channel,,  a  ralattonahip  between  die 
input  aignal-to-noiae  ratio,  and  the  multiple  coher¬ 
ence,  aimilar  to  the  one  in  tha  laat  aection  for  two 
channel,,  can  be  derived.  The  output  power  epectral 
danaity  can  again  ba  written  in  term,  of  thia  input 
aignal  epectral  deneity,  the  unknown  channel  tranafar 
function*,  and  the  affactiva  channal  noiae  aa 

S„(oi)  .  10,(01)  f  Su(o»)  4  Sv,(0»  (4. 1) 

IS^fhi)!2  •  |G,(ai)|2  fyei)!2  Su(oi)2.  id)  (4.2) 

It  ahould  ba  noted  that  (Gj(ei)|2  S  (w)  ia  tha  aignal 
power  danaity  in  tha  output  of  thcuith  channel  and 
S  .(w)  ia  tha  noiaa  powar  apactral  danaity  in  tha 
({•reel.  The  genera)  power  apcctral  deneity  func¬ 
tion  caa  than  be  written  aa 

S  (oi)  »  S  (01)  c‘(0l)  CT(oi)  4  D(W)  ,  (4.3) 

XX  u 

wbtrt 

ia  a  diagonal  matrix  with  element*  S  .(u>). 

VI 

and 

CT(ei)  ■  [C,(e>)  G^(ei) . . .  GM(ei)l  .  (4.4) 

Tha  inveree  required  to  calculate  tha  multiple  coher¬ 
ence  from  aquation  2.  2  can  now  ba  calculated  by 
using  the  following  matrix  invereion  lemma: 

{A4X*xr]'l.A*,-A-Vi4xV1x*  'xV 

(4.S) 

Uelng  thia  lemma,  tha  invaraa  of  the  apactral 
danaity  matrix  can  be  written  (aeeumtng  all  required 

invaraa,  axiat)  aa 

S"i<eo)  a  D'Soi) 

99 

-Su(olD"l(oi)  G*(oi)[l4Su(e»GT(oi)D'1(oi)G*(o>^'1 

.CT!oi)D"1(o>)  . 


(4.  b) 
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The  structure  of  this  invirn  c&a  be  men 
r1*>arly  Uy  noting  that 

T  .1  .  v  K=jt“)l**V“> 

Su(w)g‘(«)D  ‘(U)G  l“)*2.“s*iS) -  <4-7> 

j.l  vj 

Thi»  Him  ii  the  (ion  e f  til  output  .ignn).te.nei*a 
power  ntiei.  Tktnfctt,  if,  ••  is  the  tvo-cbtiuul 
case,  «•  dtfiM 

<3/N)t  .  IGjIM)!2  SjMj/S^Ul)  (4.»> 

the  inverac  of  the  bracketed  term  in  equation  4. 6  can 
be  written  aa 


^1+Su(w)GT(«)0 


-i  r  m 

■1(»)G*w>j  ■  1/  (S/Nlj  (4.9) 
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With  this,  the  ith  diagonal  eltmeat  of  S  (0)  it 
given  by 


14^  (S/N^  -  (S/N)( 

-  -<‘l _ 

M 

I  4  £  {S/NIj 

j-1 


Using  thia  and  aubatitutinf  equation  4.10  into  equation 
2.  2.  we  find 

Jr  ? 


•  1  -  1  /  Su(w)  S“(u) 


(S/N)i[  ^'(S/N)j  -<S/N),i 

- r — 33 - : - =f  ,4-ll> 

[l  ♦  (S/N)(j  jit  2,  <S/N)j  *  (S/N)j  1 

L  j*l  J 

Several  apeelal  caaea  are  of  intcreet. 

Tint  cenaid er  the  aituatien  in  which  tha 
signal *to-noi.c  ratioa  in  all  channala  are  the  aatne: 

(S/N).  •  (S/Nlj  a  (S/N)  .  (4.  12) 

Thia  give#  the  coherence  of  channel  i  with  reepeet 
to  the  other  M*1  channala  aa 

,  .2 

yi:l,2....i-l.i  +  l,...M‘ 

(S  ’S)Z  (M.l) _  ... 

*  fl  *  (S/N)]  [1  +  (M-1)(S/N)]  '  ' 


Note  that  lor  M  equal  to  two  ar  in  auction  )  the 
coherence  it  given  by 

1,1,11  ■  ■  irf^r  ■ 

However,  if  M  becomea  very  large,  the  coherence 
goea  to 

,yi!l.2.... 1-1.141,. ..M^rrMSTNiT  '  (4’,5> 


The  formal  requirement  for  thia  to  he  valid  ie  for  the 
signal.to-noisa  ratio  and  number  of  chesaela  to 
aatiofy  the  following  inequality: 

(M  -  1)  (S/N)  »  1  .  (4. 16) 

However,  haeed  on  data  from  section  3.  equation 
4.  IS  ia  identical  to  the  coherence  of  two  channala 
when  one  ha  a  an  infinite  eignal.to.noiee  ratio  mad  the 
other  haa  a  algnal>to>noiaa  ratio  (at  fraquancy  01)  of 
(S/N).  In  thia  tonae,  a  large  enough  number  of 
weak  channala  [eignal.to.noiee  ratio  of  (S/N)]  la 
equivalent  to  tha  aum  of  one  noiee.free  channel  and 
one  weak  channel. 

Tha  aecond  epeciel  caee  for  equation  4. 11  io 
whoa  tha  ith  channel  haa  a  very  large  eignal.to.noiee 
ratio.  Lotting  (S/N){  bacome  large  in  equation  4. 11 
and  keeping  all  other  eignal-to.noiae  ratioe  equal  to 
(S/N)  we  find  that 


lr  r 

7l:l,  w, . . .  i-1, 141, , , .  M 

- -  fM-I)JS/N] _ _ 

(S/N).—  *  [1  +  (M-i)  (S/N)] 

Note  that  for  low  signal -to -noiae  ratioa,  i.e. , 

(M-l)  (S/N)  «  1  . 


tha  coherence  goea  up  linearly  with  the  number  of 
channala  each  ia  considered  to  have  the  eame 
oignal-to-noiee  ratio  aa  all  others.  1.  a. ,  (S/N).  Aa 
M  becomes  larger  or  «i  ’ 

M(S/N)»  1  .  (4.19) 

thia  coherence  goea  to  one  a  a  it  would  in  tha  case  of 
two  noiaa-fraa  channala. 

Next  consider  the  caae  where  one  channel  other 
than  the  ith  channel  haa  a  very  high  signal -to -noise 
ratio  ralative  to  the  others: 


(S/N>k»>  T  fS/N)(  -  fS/N)1  -  <S/N)k  .  (4.20) 
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L'ndcr  th.ae  condition*,  equation  4.  11  it  approx) 
ninthly  *iv«n  by 


[l+(S/N),][W(S/N)k] 


,2  (S/NUS/Nl 

r  t:l ,  2, . . .  i-l,  i+1, . .  ,  Ml  pTrs7N-t’lfw,S/N)k]  <4-211 
or 

|*>:),2...|.l,l4l....M  |  '“Kilt  I  .  (4.22) 

<>  if  all  other  channel*  were  not  uaed.  If  (S/N)  la 
approximately  equal  to  (S/N). ,  thia  mean*  that  all 
weaker  channel*  could  he  neglected  and  only  the  two- 
channel  coherence  between  the  two  etronper  channele 
could  be  ueed,  Aleo  conaider  the  cae*  when  all 
channel*  including  the  ith  have  a  much  lower  aignal- 
to-noiae  ratio  than  the  kth  channel,  i.  e. , 

(S/N)k  -  -  (S/N).  .  (S/N)  all  i  /  k  (4.  25) 


Then,  while  the  coherence  of  the  ith  channel  given  the 
other*  ia  provided  by  equation  4.  21,  the  coherence 
of  the  kth  channel  given  the  other*  1* 


,2  (M-1)(S/N)  (S/N) 

|  kt],2...k-l.k+l....Ml  *(7+(s/N)kl[U(M-l)(S/N)3 


Giving  for  thia  caae 


|rk:l,2 . k-l,k*l....M  |2 


(M.ljflKS/NlI  |v  .  -  j2 

n -TnT-  1  ksVnT]  I  •  . i-i, i+i, . •  •  m|  . 


For  the  caae  of  weak  aignala  in  the  other  channel* 
(from  eouanor.  4.  22): 


The  coherence  between  the  atrong  aignal  and  the 
weaker  onea  (sea  up  linearly  with  the  number  nf 
weaker  rUennela,  '11,1a  „>-« ■„  |l,e|  «i,e  la,  peat  ,.r 

the  M  multtpl*  coherence  value*  will  be  the  on*  in 
(4.21)  the  largeat  aignal-to-noie*  ratio  channel  i* 

uaed  aa  the  reference,  which  t*  aa  expected. 


5.  A  Sample  Statiatic  for 
Multiple  Coherence 

The  true  multiple  coherence  of  a  eat  of  time 
aerie*  le  a  function  of  the  underlying  atatiatiee  of 
the**  pro******.  The  atatiatica  are  generally 
unknown  and  muet  be  estimated  from  (ample  realtea- 
tiona  of  the  procaeaea.  The  eatimatee  of  the  baetc 
atatiatica  can  then  be  ueed  to  provtd*  eatimatee  of 
the  multiple  coherence  of  the  M  underlying  etochaetic 
process**. 

The  method  of  obtaining  eatimatee  for  true 
multtpl*  coherence  it  a*  follow*.  Uaing  well  docu¬ 
mented  techniques1'  obtain  sample  estimate* 

for  each  element  of  the  croeapower  spectral  deneity 
matrix.  From  thee*  sample  estimates 

[*!,(-)  *12(W>  •"  *lM,h,»  ] 


^l*"1  S2Z,U> 


S2M<"> 


'ktl ,  2, .  • .  k-1,  k+l , . . .  M| 


(M-liriaS/N)  Jr 
[1  +(M-  I  )(S/N)]j  ,! 


(4.26) 

Further  simplify  equation  4.  26  to  the  special  case  of 
(M*l )(S/N)  <<1.  (4.27) 
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one  calculates  the  sample  estimate  for  multiple  co¬ 
herence  in  the  following  manner 


Pi:l,  2 . i-1,  i+1. . . .  M  «1  -l/(Sti<w)Sl  (u>)].  (5.  2) 

where  S^(a>)  is  the  ith  diagonal  element  of  the  inverse 
of  the  M  *by-M  sample  spectral  density  matrix  S^lw). 

Detail#  of  how  to  form  such  estimates  are  discussed 
at  length  in  the  literature.  Since  these  estimates  are 
random  variables  there  has  been  considerable  study  of 
their  distribution.  The  distributions  of  these  cross- 
and  ajitopower  spectral  estimates  are  Known  in  closed 
form ' . 
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TK«  cloa»«Uform  lor  tho  multiple* 

..  •  *v>-  -  (#  T(it«  rii|>tp«»n»< 

t:.v  uii^u  wi  v«luiti  ol  tho  multiple -coherence  toot 
fuiiitie  and  tho  relative  probability  ol  U«  bolng  in  * 
pirticuUr  bud.  All  values  art  ol  couroo  bounded  by 
zero  and  ono.  Tho  density  function  is  conditioned  on 
t no  total  number  ol  dlHerent  time  records,  or  dll* 
Ierent  stochastic  processes,  available  (p).  It  is  also 
conditioned  on  the  number  ol  independent  samples 
available  from  each  ol  tho  time  records  (N).  Thus 
the  density  function  of  the  sample  estimate  for  coher¬ 
ence  given  the  true  coherence  is  given  by* 
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Other  npniitoni  for  thia  density  valid  to  larga  values 
of  N  and  specific  rangaa  of  trua  coharanca  |  y|*  and  y 
ara  daacribad  in  reference  1.  Tharacurvaa  of  prob¬ 
ability  of  datoction  varaua  probability  of  falaa  alarm 
ara  praaantad  aa  ara  curvaa  of  probability  of  detec¬ 
tion  varaua  trua  coharanca  tor  fixed  valuaa  of  prob¬ 
ability  of  falaa  alarm. 


p-2„  ,N-p 
y  0-y)  K 


2r,(N.N.p-lIy|yJZ) 


P _ (y/N.  p,  |yP)«0  y  >  1  or  y  <0. 


In  aquation  5.  3, 


jfjf  )  ia  tha  hypargaomatric 


Thia  cxprcaaion  for  tha  danaity  function  of 
multiple  coherence  ia  both  axpanaiva  to  calculate  and  ^ 
generally  numerically  ill  Conditioned.  Thua  to  aval, 
uate  the  denaity  numerically,  additional  manipulatione 
ara  required.  Great  difficulty  can  be  encountered  in 
attempting  to  uee  computer  library  expreaaiona  for  ; 

the  hyperg cometric  function. 

For  low  valuea  of  N  and  p  we  uee  a  tranaform- 
ation  given  in  reference  12:  , 

W 

2ri(\.N.  p-li  i  yf2  y) 

*  U  *|yl2y)i>’1.'2Nj,r1(p-i-N,p-ii|r;2y).  (5.«) 

For  the  caeca  of  intereat,  (p-l-N)  ia  a  negative  7 
imager  ao  that  a  finite  aeriaa  expaneion  for  thia  latter 
nypergecmetric  function  ia  available: 

2r.(.n.N.  p-l;|y|2y) 

B 

|  |2  l-l-lN^’l-Nvp-nf  (l/yl3 

*,Hrry)  *  (p^TT  J7  ~  '»•*> 

J  ? 

where 

(p-1)  •  (p+j-lM(p-l)l  «  (P*j*i)(P*j-3).  •  .(P-1)  (*•*) 


(-N-p-1).  »  (-N»p-1  )(-N»p). . .  (-Nrp-2+j). 

U.ing  thcec  axpraaaiena  for  the  hypergeometric 
tunction  we  can  writ*  lor  y  between  taro  and  one: 
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ABSTRACT 

Thu  paper  deals  with  lire  problem  of  sound  speed  csiuna- 
uoo  as  an  aid  so  improving  our  knowledge  of  the  position  of  a 
dnflini  sonobuoy.  This  u  accomplished  by  processing  the  out¬ 
puts  of  a  spatially  displaced  set  of  hydrophone  sensors  in  an 
es tended  Kalman  filler  The  tiller's  state  vector  may  be  aug¬ 
mented  to  provide  estimates  of  the  sonobuoy i  position, 
velocity,  and  the  sound  speed  from  the  sonobuoy  to  each  of  the 
hydrophones. 

It  is  shown  that  with  jusi  timc-of-amval  difference  and 
dopplcr  difference  measurements,  system  observability  n  mar¬ 
ginal  -  directly  resulting  in  biased  estimates  of  the  sonobuoy's 
location.  With  the  additional  use  of  bearing  measurements,  n  is 
then  shown  how  this  enhances  system  observability  and,  conse¬ 
quently.  estimation  performance. 


I.  INTRODUCTION 

In  this  paper,  we  are  concerned  with  the  accurate  location 
of  a  sonobuoy  over  a  specified  interval  of  lime .  It  is  assumed  to 
be  suDiected  to  both  random  and  deterministic  forces  over  Hus 
same  tune  penod.  Because  of  this,  our  knowledge  of  the  sono¬ 
buoy's  location  is  uncertain  and  can  be  in  error  if  we  employ 
conventional  dead  recaoning  methods 

Another  approach  more  accurate  than  the  above  would  be 
to  place  a  source  emitting  a  continuous  pseudorandom  pusaian 
ugrul  aboard  the  sonobuoy.  A  set  of  three  or  more  Passive 
hydrophones,  geographically  separated  from  one  another,  can 
then  be  used  to  pick  up  the  sonobuoy  's  signal  to  provide  meas¬ 
urements  of  sound  umc-of-amva!  difierence  and  dopplet  dif¬ 
ference.  As  an  example,  one  configuration  involving  three 
hydrophones  is  shown  below  in  Figure  1. 

In  this  figure  Rj,  i»  1 , 2, 3  represents  the  distance  from  the 
sonobuoy  to  each  of  the  hydrophones  and  Cj,  i*  1 ,  2. 3  represents 
■he  sound  speed  to  each  hydrophone.  Note  that  we  are  lusng 
three  distinct  sound  speeds  instead  of  a  single  one.  This  is 
because  sound  speed,  in  a  thermally  non-homogeneous  medium, 
u  related  m  a  complex  fashion  to  such  variable!  as  time-of-day. 
distance  between  hydrophone  and  sonobuoy,  bearing  of  hydro¬ 
phone  with  respect  to  sonobuoy,  etc. 

Now  if  we  consider  any  pair  of  hydrophones,  pair  1-2  for 
instance,  then  we  can  obtain,  through  cross-correlation  methods, 
two  measurements,  time-of-emval  differences  and  doppler  differ¬ 
ence.  that  arc  related  to  the  sonobuoy's  location  by  the  following 

equations 


ft.  ft, 

ft.  ,  m  —1  .  —S. 

12  C1  c2 


(2) 


In  Equation  <  1  >.  T|.->  represent*  the  sound  timcof-amval 
difference  for  the  hydrophone  array  pair  1-2.  >1.2  in  Equa¬ 
tion  (2)  represents  the  normalized  doppler  difference.  ^ .  where 
f  is  the  base  frequency  of  the  narrowband  source. 


Using  sphencal  geometry,  the  distances  R,  and  tune  rate- 
of-changc  distances  ft,  arc  related  to  the  sonobuoy's  location  by 
the  following  equations: 


Rj  •  cos'*  | sin  *|  sinX,  ♦  cos  aj  cosX,  cos  (a^  *  8j)l  Of 
.  a3  Isin  cos  Xjoostato  -Rjl-  costj  sin  Xtl 


♦  cm  cosX,  sin  Ut  -  6,1 
sin  R, 


<41 


where 


i*l,  2.  3.  4  represents  the  sonobuoy  latitude,  longitude, 
latitude  me.  and  longitude  rate,  respectively. 

Xj.  6,  represents  the  latitude  and  longitude  of  the  i® 
hydrophone. 

In  Equation  14)  it  was  assumed  that  each  of  the  hydrophones  was 
stationary. 

Because  of  the  nonlinear  feature  of  the  above  measurement 
equations,  it  is  very  difficult  to  solve  explicitly  for  the  sono¬ 
buoy's  state  vector  (a  13,  *3,  xq).  In  the  next  section,  it  will  be 
shown  how  this  state  vector  may  Be  obtained  h  a  recursive 
manner  using  an  extended  Kalman  filter 

Following  this  we  will  show  that  when  the  sound  speeds  ere 
imprecisely  known,  then  estimates  of  the  sonobuoy's  position 
and  velocity  become  biased.  Section  3  then  presents  a  method 
for  augmenting  the  sound  speeds  to  the  sonobuoy  state  estimate 
and  it  is  diown  that  the  augmented  state  becomes  marginally 
observable.  This  directly  results  in  biased  estimates  for  the  sono¬ 
buoy’s  state  vector.  We  then  show  that  the  additional  use  of 
bearing  measurements  are  seen  to  enhance  system  observability 
to  the  point  where  estimation  of  the  sonobuoy  state  vector 
becomes  1  nbiasail. 


Finally,  in  Section  4,  we  will  summarise  our  mulls 
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Z  ESTIMATION  OF  SONOBUOY  STATE  VECTOR 


It  »  apparent  from  Equations  (I )  to  (4)  that  it  would  be 
difficult  to  solve  for  the  sonobooy  state  vector  ui  terms  of  the 
measurement  ahJ  A  recursive  approach  may  be  sought  that 
avoids  this  difficulty.  The  approach  we  use  here  involves  appli¬ 
cation  of  an  extended  Kalman  filter. 

In  order  to  implement  this  filter,  we  need  in  addition  to 
the  ineaauivmenl  Equations  ( 1  Mo  (4).  a  model  lor  slsc  dynamics 
ol  the  system.  This  is  easily  obtained  when  we  assume  a  constant 
speed  and  constant  course  drift  and,  ui  addition,  a  random  per¬ 
turbation  to  account  for  modeling  uncertainties,  environmental 
forces,  etc.  In  dtscrete^time  form,  the  equations  for  the  system 
dynamics  are  given  by: 


“l  0  At  O- 
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0  1  0  At 

0 
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0  0  10 

wj  (k) 

0  0  0  1 

w4(k) 

— 
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x:(k) 
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is  the  sonobuoy  stale  vector 


At  *  sampling  interval 

Wj  (k).  w4  (k  >  are  tero  mean,  random,  white  noise 
sequences 

E  |  w(k)wTGi|  *  Qjj  6kj  is  the  Kroneker  Delta. 


represents  the  nonlinear  dependence  of  ^  on  ilk)  defined  by 
Equation*  ( I )  through  (4)  Rk  n  the  covariance  matrix  of  an 
additive  white-noise  measurement  sequence.  H  |i(ktlAI)  is  the 
measurement  equation  that  is  linearized  about  the  most  current 
state  estimate  and  is  defined  by  : 


H(x(k4|/k» 


Bh(i(k» 

“alUkj 


xtk)  «  fitkvl/k) 


This  linearuation  was  found  to  accurately  describe  the  func¬ 
tional  characteristics  of  liie  nonlinear  measurement  model  as 
H(  • )  was  relatively  insensitive  to  variation  in  i<k). 

These  equations  were  implemented  in  the  case  where  we 
knew  the  sound  speeds  to  each  hydrophone.  Under  these  circum¬ 
stances,  the  resulting  estimation  performance  was  seen  to  be 
excellent 


However,  for  the  case  involving  unknown  sound  speeds,  the 
estimation  performance  degrades  rapidly  when  we  use  incorrect 
values  for  the  sound  speed*  in  the  Kalman  filter.  An  example  of 
this  is  shown  in  Figures  2  and  3  for  a  case  involving  an  average 
sonobuoy  drill  velocity  of  1  knot  and  course  of  4  3  degrees.  In 
Figure  2,  the  solid  curves  represent  motion  of  the  truth  mode! 
and  the  dotted  curves  represent  estimates  of  the  sonobuoy  state 
as  provided  by  a  Kalman  filter.  The  measurements  were  assumed 
noisy  with  a  standard  deviation  in  t  and  a  of .  1  second  and 
1  x  10*5.  respectively.  The  estimation  error  and  a  priori/a  pos¬ 
teriori  standard  deviation*  are  presented  in  Figure  3.  A  set  of 
three  hydrophones  were  configured  essentially  as  vertices  of  an 
equilateral  mangle  with  the  sonobuoy  starting  near  the  center  of 
this  triangle.  The  truth  model  sound  speed*  were  set  at:  C|  * 
4857  ft/sec,  C2  “  4957  ft/scc,  and  C3  *  4757  ft/sec  around 
nominal  values  from  Unck  ( I  ].  The  filler  used  an  incomct  value 
of  5000  ft/sec  for  all  three  speeds  and  no  attempt  to  estimate 
them  was  made.  It  is  clear  from  Figure  2  that  significant  biases 
in  estimated  position  and  speed  develop. 

In  the  next  section,  we  will  see  ho*  the  use  of  bearing 
measurements  enhances  system  observability  and  allows  unbiased 
estimation  of  the  sonobuoy  stale  vector. 


If  we  linearize  the  measurement  equation  about  our  most  recent 
estimate  of  the  state,  &(k),  to  obtain  a  linear  measurement 
matrix,  H(k(k)),  we  can  then  implement  the  extended  Kalman 
filter  ui  real  time.  In  vecloc-matnx  form,  these  equations  are 
recursively  given  by  (3) ; 

Mk+l/k)  »  ♦(At)  t<k/k)  (6) 

Plk+l/k)  *  AlAIIPlk/k)  ♦T(AI)  +  Qk  (7) 

Hk+l/kel)  *  l(kel/k|eKkUk-h(*<k»t/k))|  («) 


3.  SOUND  SPEED  ESTIMATION 

Because  of  the  bum  that  develop  when  we  use  inconeci 
round  ipeedi  in  the  Kalman  Tiller,  n  was  fell  that  estimation  of 
these  round  speeds  in  addition  to  the  sonobuoy  states  would 
eJiminaic  the  biases.  This  is  done  amply  by  augmenting  the 
sound  speeds  to  the  original  four  element  sonobuoy  suie  vector 
The  assumption  wet  made  that  the  sound  speed  was  constant  and 
unknown  over  a  sufTicienily  long  time  interval.  Its  SUte  equation 
b  then  given  by: 

«i  *  0;  i  •  1.2.  J  (II) 


Plk+l/kel  I  -  |l- Kk  H(i(k*l/k)l)  Ptt+l/k)  (9) 


Kk  •  Plkvl/k)  H(l(k*l/k)>  |H(8(k*l/k))  Plkvl/k) 
HT(Jlk»l/k))»  Rkl  1 


(10) 


Equations  16)  and  (7)  represent  predictions  of  the  stale 
estimate  and  covanance  lo  tune  Ik-1),  whereat  bquationt  (I) 
through  (I0|  update  the  state  estimate  and  covanance  to  account 
for  a  new  measurement.  In  the  Jailer  act  of  equations,  the 
measured  sector  is  given  by  ak '  •  I T  ^(kl,  a^k)! .  where  h(  • ) 


Implementation  of  the  extended  Kalman  filler  for  the  aug¬ 
mented  state  vector  is  therefore  rathet  simple  and  straightfor¬ 
ward  and  will  not  be  presented  here. 

Upon  implementing  the  augmented  stale  filter,  it  was  found 
that  biases  still  developed  in  I  he  sonobuoy  stales  and  round 
speed!  for  a  number  of  tcenanos  involving  different  surfing 
positions,  speeds,  courses,  etc.  For  this  reason,  the  observability 
of  Ihc  linearized,  augmented  tlate  system  wai  explored  using  the 
information  matrix  |2|.  In  Ihc  tilustion  involving  no  process 
none  and  Hale  vector  a  pnon  information,  the  information 
malna  ■  equal  to  the  inverse  of  Ihc  Kalman  filler  COVS nance 
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■Minn,  H  (k/k  i  Thu  uuiru  imiil  be  pouim  dcTuult  for 
u<KlM>iu  obwnikkiy  and,  providad  tht  abort  caoditaoM 
apply .  u gate*  an  rcoanan  (dan  bp: 

r'(k/k)»  •Tl-Al>fl  (kl/k-l)  ♦(-£!* 

<I2> 

♦  Mt  (I|k*l»  Rk-1  HHIk-IMir1  (0/0)  >0 

•here  ♦(  At)  is  the  stale  transition  malru  defined  in  (5), 

Hi  It  a- J )» n  the  measurement  mains  Unearned  about  the  stale 
vector  S(M). 

To  aiaeia  the  property  of  tlochasiic  observability,  the 
normalized  eigenvalues  of  this  matrix  were  computed  I  nor* 
maliicd  to  one)  and  plotted  as  a  function  of  lime.  Figure  4  shows 
the  results  that  were  obtained  for  the  caae  involving  one 
unknown  and  constant  sound  speed  One  of  the  sonobuoy  posi¬ 
tion  eigenvalues  becomes  ill-conditioned  and  exhibits  a  smaller 
maximum  magnitude  than  the  other  position  eigenvalue  by  a  few 
orders  of  magnitude.  Thu  analysis  was  repeated  for  e  number 
of  different  sonobuoy  starting  positions  and  headings,  and  the 
same  genera)  result  was  obtained,  ill-conditioned  behavior 
of  one  of  the  eigenvalues.  Because  of  this  and  the  fact  that  the 
state  oumates  were  biased,  we  concluded  that  the  system  was 
marginally  observable 

As  a  means  of  improving  the  property  of  system  observa¬ 
bility.  we  introduced  an  additional  measurement  of  bearing  to 
use  m  the  Kalman  filter.  We  assumed  that  it  was  contaminated  by 
additive  white  muse  with  a  standard  deviation  of  .1  degree.  These 
measu/emenU  defined  the  bearing  of  the  tonobuoy  wtth  reipcct 
to  each  hydrophone,  and  are  normally  available  as  outputs  front 
the  same  set  of  hydrophones  without  introducing  additional 
sensors 

The  Kalman  filter  was  then  rerun  with  these  additional 
measurements  and  the  sonobuoy  tracking  performance  was 
found  to  be  very  food  Figures  5.  b,  7.  and  b  summarize  the 
result*  of  one  case  where  wc  estimate  the  four  sonobuoy  states 
and  three  sound  speeds  In  Figure  5,  convergence  Of  the  sono¬ 
buoy  state  estimates  to  their  true  values  is  seen  to  be  very  rapid. 
Noise  was  introduced  in  the  sound  speed  truth  mode)  to  add 
more  realism  to  the  problem  and  this  variation  u  reflected  by  the 
variation  in  sound  speed  of  the  truth  model  as  shown  in  Figure  6. 
Figures  7  and  6  present  estimation  error  and  standard  deviations 
of  the  sonobuoy  states  and  sound  speeds.  These  errors  are  seen 
to  fall  well  within  the  t  two-sigma  value  as  one  would  expect 
for  unbiased  estimation  A  number  of  other  cases  were  also  run 
and  arc  not  presented  here  because  ol  space  Imitations.  For 
l he ir  oilier  cases  it  was  found  that  tracking  performance  was  as 
guud  as  belore  and  all  of  the  earlier  noted  biases  in  the  state 
estimates,  when  processing  only  time-delay  and  doppter  differ¬ 
ence  measurements,  had  disappeared. 


4.  CONCLUSIONS 

To  summarize  the  above  results,  we  saw  that  to  estimate  the 
sonobuoy  position  and  velocity  accurately  wc  need  to  also  asti- 
matc  the  sound  speeds  to  each  hydrophone  when  they  art 
unknown.  With  the  aid  of  the  information  matrix  we  found  that 
the  Unearned  system  was  margmaJiy  observable  when  attempting 
to  estimate  both  the  sonobuoy  states  and  sound  speeds. 

Because  the  linearized  measurement  model  accurately  tepee* 
tented  its  nonlinear  counterpart,  we  concluded  that  the  nonlinear 
system  was  also  marginally  observable.  This  eras  also  reflected 
by  the  biased  estimates  of  sonobuoy  position  when  implementing 
the  Kalman  filler. 

Finally,  by  incorporating  and  additional  independent  mea¬ 
surement  (bearing)  the  system  observability  was  enhanced  con* 
siderably  and  we  weir  then  able  to  track  the  sonobuoy 's  motion 
quite  accurately. 

Although  not  reported  in  the  last  section,  wc  also  looked 
at  the  use  of  time-delay  and  doppler  difference  measurement* 
of  more  than  three  hydrophones  tas  opposed  to  time  delay, 
doppter  difference,  and  bearing  measurements  of  just  three 
hydrophones)  as  a  means  of  improving  system  observability  In 
general,  it  was  found  that  if  one  had  n  hydrophones  with  n 
unknown  sound  speeds  from  the  sonobuoy  to  each  hydrophone, 
then  wc  could  at  most  estimate  the  sonobuoy  stale  vector  and 
(n*l )  sound  speeds;  one  of  these  sound  speeds  had  to  be  known 
exactly  when  implementing  the  Kalman  filter. 

At  tht  preaent  time,  we  arc  continuing  research  in  this  area 
of  aonobuoy  state  and  sound  speed  estimation  to  more  fully 
explore  the  sensitivity  of  the  estimation  process  to  geometrical 
placeman!  of  the  hydrophones,  lime  between  measurements,  etc. 
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Hydrophone  1 


Figure  1.  Hydrophone /eonobuoy  geometry. 
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CHNDOT  -  APRIOR  STD.  DEV. 

Figure  3.  Sonobuoy  estimation  error  and  standard 
deviation— incorrect  sound  speeds. 
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••LONG.  RA1  E  -  2.767+01 

SOUND  SPEED  *  1.420+00 
COURSE  -  90 

Figure  4.  Normalized  eigenvalues  of  the  information 
matrix. 
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Figure  5.  Estimation  of  sonobuoy  states. 
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Figure  8.  Sound  speed  estimation  and  standard 
deviation. 
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OVERVIEW  OF  INTERSENSOR  CORRELATION  FOR 
OCEAN  SURVEILLANCE _ 


In  a  typical  large-area  ocean  surveillance  situation,  data  is 
generated  by  many  different  sensors  due  to  the  presence  of  various 
surface  ships,  aircraft,  etc.  If  the  data  collected  by  the  sensors  is 
overlaid  in  a  common  coordinate  system,  then  a  picture  such  as  that 
shown  in  Figure  1-1  below  results.  The  under  lying  assumption  of  the 
picture  is  that  each  individual  sensor  has  formed  a  picture  of  the  sur¬ 
veillance  area  based  only  on  its  own  data.  In  some  cases,  this  may 
require  the  solution  of  a  multitarget  tracking  problem  [l]  at  the 
sensor  level.  In  Figure  1-1,  for  example,  sensor  S2  must  discrim¬ 
inate  and  track  the  two  closely  spaced  targets  T1  and  T2. 

When  the  multitarget  tracking  problem  is  solved  separately 
for  each  individual  sensor,  a  somewhat  redundant  view  of  the  sur¬ 
veillance  area  may  result.  The  different  data  collection  systems 
SI,  S2,  ....  may  report  on  the  same  target.  When  this  situation 
arises,  the  final  step  in  the  process  of  forming  a  complete  picture 
of  the  surveillance  area  is  to  perform  inter  sensor  correlation. 
Referring  again  to  Figure  1-1,  we  see  that  target  T1  is  seen  by 
sensors  SI  and  S2.  The  tracks  of  T1  produced  by  the  two  sensor 

'  *  t 

systems  will  not  overlay  exactly- -thus  the  requirement  for  some 
sort  of  decision  process  (automatic  or  manual)  to  accomplish  inter¬ 
sensor  correlation. 

This  paper  is  devoted  exclusively  to  automatic  (i.  e. ,  com¬ 
puter)  algorithms  for  carrying  out  the  tasks  associated  with  inter¬ 
sensor  correlation.  One  particular  algorithm  (the  "k-track  cluster¬ 
ing"  algorithm)  is  discussed  in  some  detail  in  order  to  highlight  the 
fundamental  combinatorial  structure  of  large  area  surveillance. 
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The  working  definition  of  Multisource  Interaction  (MSI) 
adopted  in  this  paper  is  one  which  is  effectively  decoupled  from  the 
problems  of  sensor  allocation  and  sensor-level  multitarget  tracking. 
By  that  we  mean  that  the  sensors'  fields  of  view  are  controlled  by 
a  separate  process,  and  that  the  multitarget  tracking  problem  is 
solved  on  a  sensor-by-sensor  basis.  The  responsibility  of  MSI 
data  processing  is  then  to  carry  out  inter  sensor  correlation  for 
various  track  file  data  bases.  The  MSI  issues  addressed  herein 
arise  due  to  the  partitioning  of  a  large  surveillance  data  base  into 
smaller,  more  manageable  segments.  As  Figure  1-2  indicates,  data 
is  generated  by  many  sensors  due  to  objects  present  in  several  spa¬ 
tial  sectors  during  some  specified  time  interval.  If  the  multitarget 
tracking  problem  is  solved  separately  for  each  individual  sensor  (see 
[  1  ] ),  then  track  estimates  for  objects  in  the  surveillance  area  become 
consolidated  first. at  the  level  of  individual  sensors.  This  is  indicated 
schematically  in  Figure  1-2  as  the  consolidation  of  the  raw  sensor 
measurements  (data  bases  A,  B,  C)  into  a  "sensor-level"  data  base 
D  that  consists  of  a  complete  track  file  for  that  specific  sensor. 

The  final  step  in  the  process  of  forming  a  complete  picture 
of  the  objects  present  in  the  surveillance  area  is  the  solution  to  the 
intersensor  correlation  problem.  At  this  level  of  the  MSI  data  pro- 
cessing  hierarchy,  the  track  files  D,  E,  F  belonging  to  the  different 
sensor  systems  are  compared  to  determine  the  correlation  between 
tracks  in  separate  data  bases.  As  indicated  in  Figure  1-2,  the  out¬ 
put  of  an  intersensor  correlation  algorithm  is  an  entry  (or  entries) 
into  a  master  data  base  G  that  consists  of  track  files  on  an  object- 
by-object  basis  for  the  entire  surveillance  volume.  If  the  intersensor 
correlation  problem  is  solved  correctly,  the  master  data  base  will 
contain  a  complete  list  of  object  tracks,  without  duplication,  and  in 
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Figure  1-2.  Overview  of  MSI  Hierarchical  Information  Processing. 
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in  the  algorithm  described  in  Section  2.  is  latitude,  latitude  rate 
longitude,  longitude  rate.  ) 

1.  1  Description  of  the  Intersensor  Correlation  Problem 


The  multitarget  tracking  data  processing  done  at  the  sensor 
level,  using  algorithms  such  as  those  discussed  in  [l],  produces 
track  files  consisting  of  object  state  vectors  and  (sometimes)  estima¬ 
tion  error  covariance  at  specific  points  in  time. 

Typically,  the  sensor  level  track  files  will  contain  information 
such  as  that  shown  below  in  Table  1-1. 


Table  1-1.  Typical  information  content 
of  an  MSI  Data  Base. 


\Data 

State  Vector  Error 

Object 

TracX 

State  Vector 

Covariance 

Name 

1 

3*1  .  3*1 

x  (tj),  x  (t2) 

v,v.v,v 

QUEEN 

MARY 

2 

8*2/+  x 

x  (t3) 

Unknown 

Unknown 

3 

6*3  6*3 

x  (t4),  x  (tg). 

6p3(t4),6P3(t5), 

Unknown 

6*3,,.  . 
x  <t6) 

V<v 

V* 


—  | 


The  notation  used  in  Table  1-1  is  defined  as  follows: 


<V 


~  sensor  level  track  file  estimate  of  the  state 
at  time  t^  of  a  target  seen  by  sensor  j.  The 
individual  tracks  are  denoted  by  the  index  i. 


^Pl  (t^)  —  the  sensor  level  track  file  error  covariance 
matrix  of  the  measurement  ^x*  (t,  ). 


It  is  often  the  case  that  the  (x,P)  data  recorded  in  the  data 
base  is  itself  generated  by  a  tracking  filter  of  some  kind.  Referring 
again  to  Figure  1-2,  recall  that  the  (x,  P)  records  are  associated 
with  sensor  level  track  files  <D,  E,  F).  Thus,  (x,  P)  records  are 
created  by  processing  the  raw  sensor  data  for  one  sensor  with  a 
filter  to  produce  a  chain  of  target  track  points. 

Because  of  the  diversity  of  sensors  that  are  operating  in  a  large 
area  surveillance  system,  the  time  points  t^  at  which  the  object  states 
(latitude,  longitude,  etc.  )  are  estimated  may  not  coincide.  Further¬ 
more,  the  number  of  time  points  per  object  and  even  the  state  vector 
coordinate  system  where  the  sensor -level  multitarget  tracking 
problem  is  resolved  will  differ  for  the  various  sensors'  track  files. 

The  diverse  nature  of  the  sensor -level  track  files  introduces  an 
added  degree  of  complexity  into  the  intersensor  correlation  problem. 

Based  on  these  last  comments,  a  more  accurate  definition 
of  intersensor  correlation  might  be  as  follows:  to  cluster  together 
those  tracks  corresponding  to  the  same  object,  and  to  produce  a 
composite  estimate  of  each  object's  motion  in  a  single  coordinate 
system  common  to  all  sensors. 
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1.  2  Algorithms  for  Intersensor  Correlation 

There  are  two  basic  components  of  any  intersensor  correla¬ 
tion  process:  a  similarity  measure  to  quantify  the  "distance" 
between  tracks,  and  a  clustering  criterion  to  carry  out  the  actual 
correlation  of  tracks.  Since  the  tracks  to  be  clustered  together 
are  in  fact  samples  from  a  stochastic  process,  similarity  measures 
can  be  based  theoretically  upon  the  distance  between  stochastic  pro¬ 
cesses.  In  essence,  we  simply  process  the  data  from  different 
track  files  with  a  single  Kalman  filter;  the  residuals  from  the  filter 
are  then  used  to  measure  the  "distance"  between  the  tracks. 

The  clustering  criteria  used  to  determine  the  intersensor 
correlations  can  range  from  various  forms  of  pairwise  correlators 
to  the  more  sophisticated  k-wise  correlation  to  be  discussed  below. 

As  an  example  of  a  simple  pairwise  clustering  criterion,  consider 
the  following: 

Step  1.  Compute  the  similarity  S..  between  all  pairs 
of  tracks 

Step  2.  If  S.j  is  greater  than  a  preset  threshold  6,  then 
Track  i  and  Track  j  are  declared  to  be  the  same 
object  (i  —  j) 

Step  3.  If  i  —  j,  and  j  —  k,  then  i  **  k  (i.e.  ,  all  three 
tracks  are  the  same  ship). 

A  little  reflection  upon  the  algorithm  presented  above  will  indicate 
to  the  reader  that  an  iterative  application  of  Step  3  will  produce  a 
complete  picture  of  the  correlations  in  any  particular  data  base.  It 
is  also  clear  that  the  algorithm  is  vulnerable  to  at  least  one  type  of 
cascading  error.  This  is  shown  below  in  Figure  1-3,  which  contains 
four  targets  (Tl,  T2,  T3,  T4)  on  parallel  courses.  If  S12'  S23’  S34 
are  all  greater  than  6,  then  the  pairwise  algorithm  outlined  above  will 
declare  all  four  sensor-level  track  files  to  be  the  same  ship. 
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Figure  1-3.  Cascading  errors  in  a  simple  pairwise  intersensor 
correlation  process. 


The  object  of  these  last  remarks  is  to  illustrate  that  simple 
automatic  algorithms  for  intersensor  correlation  may  have  hidden 
difficulties,  and  that  there  is  some  merit  in  working  with  the  more 
complex  decision  processes  discussed  in  the  following  sections. 

1.2.1  Track  File  Similarity  Measure 

This  section  will  describe  a  simple  technique  for  measuring 
the  similarity  between  different  tracks  in  an  MSI  data  base.  For 
example,  in  Table  1-1  above,  a  fundamental  question  that  the  analyst 
might  wish  to  ask  is  the  following:  how  similar  is  track  #1  (QUEEN 
MARY)  to  track  #3?  The  basic  technique  we  employ  is  illustrated  in 
Figure  1-4  below.  First,  the  points  t^,  t^,  t^,  t^,  t^,  t^  are  ordered 
in  terms  of  increasing  time.  Following  this,  the  sensor -level  state 
vectors  are  input  to  a  Kalman  filter,  and  the  residual  sequence  is 
monitored. 


3a1  .  6a3 

x  (t  ),  x  (t2),  x  (t4 

6  a  3  .  6  a  3 

x  (t5),  x  (tfe) 


Figure  1-4.  Measuring  the  similarity  of 
different  track  files. 

When  testing  the  similarity  of  two  tracks  i  and  j,  it  is  neces¬ 
sary  that  the  test  be  carried  out  in  a  coordinate  system  common  to 
both  tracks.  Since  the  tracks  i  and  j  come  from  different  sensor- 
level  track  files,  it  is  possible  that  they  are  represented  as  vectors 
x  of  different  form.  Asa  simple  example,  the  state  of  track  i  might 
be  represented  as  latitude /longitude  only,  while  track  j  might  be 
represented  as  latitude,  longitude,  course,  and  speed. 
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In  the  algorithm  described  in  this  paper,  sensor-to- sensor 
track  correlation  is  carried  out  for  all  tracks  in  a  single,  common 
coordinate  system.  We  assume  that  target  motion  is  adequately  modeled 
by  a  stochastic  difference  equation  of  the  form 

y{k+l)  =  ?(k+l,  k)  y{k)  +  v(k),  k  =  0,  1,  .  .  . ,  (1) 

where 

y(k)  ~  vector  describing  object  state  at  time  t^  (latitude, 
latitude  rate,  longitude,  longitude  rate) 

tp{- ,  • )  ~  state  transition  matrix  model  (for  "straight-line" 
motion) 

v(k)  -  noise  model  compensating  for  inaccuracies  in  mod¬ 
eled  track  motion  (e.g.,  small  random  course  changes) 

The  state  vectors  x  from  the  track  files  are  treated  as  if 
they  are  "measurements"  of  actual  target  motion;  the  error  covar¬ 
iance  matrices  P  of  the  state  vectors  are  treated  as  though  they 
are  "noise"  in  the  measurements.  Using  this  pseudo-measurement, 
we  have; 

jx\k)  =  cj(y(k),k)+V(k)  (2) 

where 

y(k)  -  actual  target  state  at  time  t^ 

cV»  •)  ~  observation  equation  for  sensor 

^wX(k)  ~  "measurement"  noise  with  covariance  ^PX(t^) 

10 


I 


Using  the  x  as  inputs,  a  Kalman  filter  can  be  used  to  esti¬ 
mate  the  true  ship  motion  y(k)  under  the  assumption  that  the  inputs 
x  from  the  various  sensor-level  track  files  are  due  to  the  same  ship. 
As  the  Kalman  filter  operates,  residuals  r(k)  are  produced  whose 
magnitude  can  be  predicted,  using  the  covariances  P  from  the  track 
files  and  the  covariances  of  v(k)  in  Equation  (1)  above.  In  particular, 
the  covariance  matrix  A(k)  of  the  residuals  can  be  computed. 

If  the  observed  variance  of  the  residuals  r(k)  compares 
favorably  with  their  predicted  value  A(k),  then  the  sensor- level 
tracks  whose  state  vectors  were  used  to  compute  the  residual 
sequence  {r(k)}  are  said  to  be  "similar.  " 

Note  that  in  the  two-track  case,  if  the  vectors  for  the  two 
tracks  i  and  j  being  tested  for  correlation  are  in  the  same  coordi¬ 
nate  system  and  at  the  same  time  point,  then  we  could  alternatively 
use  euclidean  distance  as  a  similarity  measure,  for  example 

2 

distance  (i,j)  =  II  x*(t)  -  *xJ(t)ll  (3) 

S.j  =  exp  (-distance  (i,  j) )  (4) 

However,  the  tracks  we  deal  with  from  the  various  MSI  files  are 
often  not  in  the  same  coordinate  system  and  not  at  the  same  time 
point.  Therefore,  we  resort  to  the  Kalman  filtering  technique 
described  above  to  compute  distance  between  tracks. 

Each  potential  correlation  w  is  the  result  of  a  hypothesis 
test  that  uses  a  likelihood  function  p(w)  determined  from  the  motion 
model  Equation  (1).  For  example,  if  Track  #1  and  Track  #3  from 
Table  1-1  are  to  be  tested  to  determine  their  similarity,  then 

J  3 a  1  3a  1 ,  6a3  6a  3  .  6a 3..  .  \ 

«  =  I  X  (tj),  X  (t2),  x  (t4),  x  (t5),  x  (tfc)J 
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is  a  2 -track  correlation  and  the  likelihood  function  is  found  by  com¬ 
puting  a  sequence  of  estimates  y(t.)  of  ship  position  using  Equations 
(1 )  and  (2)  and  «.  The  rule  that  decides  that  tracks  i  ,  i  i  ... 

X  Ct 

are  correlated  is  then 

S.  ,  i  =  In  p  (w  [  f  y  (t  ) }  )  2  6.  (5) 

Based  on  the  log  likelihood  decision  function,  the  set  of  potential 


correlations  is 


D  =  |  u  I  in  p(a>|{y{tk)}  ^  . 


As  noted  above,  the  actual  computation  of  S^^ij...  *s  car* 
ried  out  by  a  Kalman  filter,  using  the  innovations  sequence 


r(k)  =  V(k)  -  cj(*(k,  k-l)y  (k-1),  k)  . 


The  negative  log -likelihood  function  is  given  in  this  case  by 


-  In  p(u)  =  y  ^  dim^xl(tk))fn  2ir  +  y  ^  fn  |A(k)| 


+  y  r  (k)*  A  (k)"  1  r(k). 


where  A(k)  is  the  covariance  of  r(k)  computed  by  the  Kalman  filter 
and  dim(^xl(t,  ))  is  the  dimension  of  Jxl(t,  ). 


The  summation  Zr(k)*A(k)  *r(k)  is  a  chi- squared  random 
variable,  as  are  its  individual  terms  r(k)*A(k)  *r(k).  The  remain¬ 
ing  terms  in  Equation  (8)  are  deterministic.  As  a  result,  the  similar¬ 
ity  test  Equation  (5)  can  be  written  in  the  form  of  the  following  cumu¬ 
lative  chi-square  test: 

n 

£  rfk/ilkf’rlk)  <  X*.  (9) 

k=  1 

2 

Note  that  X  is  a  function  of  the  number  of  individual  measurements 
c 

in  the  correlation. 

The  actual  construction  of  the  complete  set  of  potential  corre¬ 
lations  D  is  accomplished  through  the  use  of  a  depth  first  back¬ 
tracking  algorithm.  By  this  we  mean  that  as  many  tracks  as  possible 

are  added  to  a  potential  correlation  before  backtracking  occurs. 

When  backtracking  occurs,  all  points  associated  with  the  track 
being  stripped  out  of  the  correlation  are  removed. 

1.2.2  Discriminating  Between  Correct  and  Incorrect 
Correlations _ 

Because  of  the  statistical  nature  of  the  decision  process  (9) 
used  to  construct  the  potential  correlation  file  D,  some  tracks  may 
be  included  in  more  than  one  correlation.  It  may  be  difficult  to 
decide  which  correlation  is  correct  (if  any)  based  only  on  the  sim¬ 
ilarity  measure.  The  Bayesian  decision  process  is  structured  to 
alleviate  this  difficulty  by  evaluating  only  complete  pictures  of  the  sur¬ 
veillance  area.  The  potential  correlations  in  D  are  matched  together 
in  every  way  possible  until  the  best  possible  "global"  picture  of  the 
area  is  found.  The  basic  constraint  operable  at  each  stage  in  this 
combinatorial  decision  problem  is  that  a  track  can  be  used  in  at 
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most  one  correlation.  The  actual  optimization  problem  is  the 
following: 
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maximize 

subject  to  the  constraint 
that  tracks  are  only  used 
once 


The  mathematical  form  of  the  problem  is  as  follows: 

,t 

max  d  v 


subject  to 


Bp  <  1 


(10) 


(11) 


(12) 


v  binary  (13) 

where  d  is  a  vector  of  similarity  measures  •  •  •  »  v  is  a  binary 

vector  aenoting  which  elements  of  D  were  selected,  B  is  a  matrix  of 
zeros  and  ones  chosen  to  enforce  the  constraint  mentioned  in  (10),  and 
1  is  a  vector  of  ones. 

t 

1.2.3  Integer  Programming  Methods 

Problem  (11)  -  (13)  is  a  0-1  integer  program  that  can  be 
solved  using  implicit  enumeration  techniques.  The  basic  idea  of 
these  techniques  is  to  minimize  through  the  use  of  appropriate  tests 
the  extent  to  which  the  Bayesian  decision  tree  must  be  examined 
before  a  solution  of  (11)-  (13)  is  found.  The  actual  technique 

employed  in  the  algorithm  is  discussed  in  some  detail  in  [  1  ] .  Since  the 
combinatorial  structure  of  the  Bayesian  decision  process  is  the  same 
for  both  k-track  correlation  and  the  multitarget  tracker  discussed  in 
[l] ,  the  algorithm  can  be  used  without  modification. 
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2.  0  NUMERICAL  RESULTS 

This  chapter  will  discuss  some  of  the  characteristics  of  the 
k-track  correlation  process  as  it  is  actually  applied  to  data.  Table 
2-1  below  summarizes  a  set  of  two  test  cases  that  were  analyzed. 
Runs  #1  and  2  were  synthetic  data  cases  in  which  intermittent,  noisy 
track  reports  were  made  on  six  ships.  The  synthetic  data  test  cases 
will  be  discussed  in  some  detail  to  illustrate  the  basic  structure  of 
the  MSI  correlation  problem. 


Table  2-1.  Summary  of  numerical  results 

using  k-track  correlation  algorithm. 


Parameter 

File 

Number  of 
Potential 
Cor  relations 

PAR  1  A 

11 

PAR  IB 

26 

Ttmt  Required  Time  Require 

for  Potential  Number  of  by  Integer 

Correlation*  Correlation*  Program 

face.  )  Selected  (sec.  ) 


Number  of 
Individual 

T rack  Reports 
in  Data  Base 

Number 
of  Data 
Point* 

10 

30 

10  1 

30 

Description  of  Simulated  Data  Bases 


The  truth  model  for  Runs  #1  and  2  is  contained  in  Figure  2-1 
and  Tables  2-2  and  2-3  below.  The  actual  data  bases’(ORTlA  and 
ORT1B)  used  in  Runs  #1  and  2  are  contained  in  Tables  2-4  and  2-5. 

The  situation  modeled  in  data  bases  ORT1A  and  ORT1B  is 
one  in  which  six  simulated  ships  are  seen  and  reported  by  a  fictitious 
suite  of  MSI  sensors  on  different  occasions  within  a  60-hour  time 
period.  The  ships  are  assumed  to  hold  constant  latitude  rate/ 
longitude  rate  courses  for  the  entire  60  hours.  The  sensors  that 
are  brought  to  bear  on  the  situation  are  of  two  types,  as  indicated 
in  Table  2-6.  Type  l  sensors  measure  latitude,  longitude,  and  their 
time  derivatives.  Type  2  sensors  measure  latitude  and  loncitude 


> 


Six- shi 


Table  2-2.  Summary  of  6-ship  simulation. 
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Table  2-4.  Synthetic  MSI  data  base  for  Run  #1 
(data  base  ORT1A). 


x - MOVEMENT  INFORMATION - X 


NORTH 

LATITUDE 

EAST 

LONGITUDE 

TIME 

LAT 

RATE 

LONG 

RATE 

1 

43200.0 

4.1492 

0 . 9206E-05 

152.9049 

-0 . 4658E-04 

r> 

129600.0 

5.0456 

0 . 8344E-05 

149.2870 

-0 . 4341E-04 

3 

43200.0 

4.2937 

0,591 IE-05 

152.9590 

-0 . 4638E-04 

4 

86400.0 

4.8076 

0 ♦ 8067E-05 

151.0732 

-0.4574E-04 

5 

172800.0 

5.5932 

0 . 9221E-05 

146.9504 

-0 . 4462E-04 

6 

216000.0 

5.8718 

0 . 5875E-05 

145.1482 

-0 . 4542E-04 

7 

86400.0 

5.5847 

0 » 0000E+00 

156.9793 

0.  OOOOE+OO 

8 

129600.0 

6.5455 

0 • OOOOE+OO 

159.5486 

0. OOOOE+OO 

9 

172800.0 

6.9731 

0 . OOOOE+OO 

162.6776 

0. OOOOE+OO 

10 

0.0 

3.8903 

0. 1803E-04 

150.9037 

0 ♦ 6873E-04 

11 

86400.0 

5.4676 

0. 1771E-04 

156.6906 

0 . 6771E-04 

12 

129600.0 

6.1707 

0. 1763E-04 

159.5567 

0 . 6677E-04 

13 

172800.0 

7.2429 

0. 1972E-04 

162.5593 

0 . 6782E-04 

14 

86400.0 

5.3501 

0 . 1733E-04 

156.6992 

0 . 6804E-04 

15 

129600.0 

6.2073 

0.1607E-04 

159.6627 

0 . 6680E— 04 

16 

172800.0 

6.9876 

0, 1758E-04 

162.7621 

0 , 6848E-04 

1  7 

216000.0 

7.7721 

0. 1728E-04 

165.5123 

0 . 6816E-04 

13 

43200.0 

1.8336 

0. OOOOE+OO 

153,6828 

0.  OOOOE+OO 

19 

172800.0 

1.8522 

0. OOOOE+OO 

149.4227 

0. OOOOE+OO 

20 

0.0 

6.1539 

0 . 2936E-04 

158.1071 

-0 . 6317E-06 

21 

172800.0 

10.6991 

0 . 2757E-04 

157.9753 

0. 1652E-05 

2° 

0.0 

6.0474 

0 . 2727E-04 

153.0660 

-0 . 4371E-06 

23 

43200.0 

7.2062 

0 . 2697E-04 

158.1080 

-0 . 3568E-05 

24 

86400.0 

8. 3819 

0 . 2803E-04 

158.1236 

0 . 6289E-06 

*r 

172800.0 

10.8695 

0.2667E-04 

158.0033’ 

-0.91 70E-06 

86400.0 

5.9870 

0 • 5672E-04 

151.9055 

-0. 1170E-04 

"7 

/ 

216000.0 

13.7665 

0 • 6102E-04 

150,7840 

-0 • 1249E-04 

28 

43200.0 

5.8855 

-0 . 7186E-04 

155.1613 

0 . 2583E-04 

"f  L. 

129600,0 

-0.6172 

-0 . 7366E-04 

157.4634 

0 . 2694E-04 

7  ( > 

216000.0 

—  6  *  9036 

-0.7588E-04 

159.7972 

0 . 2731E-04 

I 


I 
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Table  2-5.  Synthetic  data  base  for  Sun  #2 
(data  base  ORT1B). 


x - MOVEMENT  INFORMATION - X 


NORTH 

LATITUDE 

EAST 

LONGITUDE 

TIME 

LAT 

RATE 

LONG 

RATE 

1 

43200.0 

2.3662 

0 . 1931E-04 

151.8230 

-0 . 5395E-04 

2 

129600.0 

5.0791 

0 . 1099E-04 

151.2165 

-0.2292E-04 

3 

43200.0 

3.8109 

-0 . 1254E-04 

152.3634 

-0 . 5202E-04 

4 

86400.0 

5.8242 

0 . 8303E-05 

151.2872 

-0 . 4569E-04 

5 

172800.0 

7.4293 

0.1 946E-04 

145.6509 

-0.3447E-04 

6 

216000.0 

7,0897 

-0. 1288E-04 

145.4386 

-0.4188E-04 

7 

86400.0 

5.8711 

0 . OOOOE+OO 

158.4497 

O.OOOOE+OO 

8 

129600.0 

8.4906 

O.OOOOE+OO 

157.9215 

O.OOOOE+OO 

9 

172800.0 

5.7781 

0. OOOOE+OO 

162.9507 

O.OOOOE+OO 

10 

0.0 

2,9031 

0.1 850E-04 

150.0369 

0.8164E-04 

11 

S6400.0 

4.6994 

0 . 1543E-04 

155.5623 

0.7047E-04 

12 

129600.0 

4.7428 

0 . 1469E-04 

158.0031 

0 . 6056E-04 

13 

172800.0 

8.4769 

0.3482E-04 

161.768 7 

0.6976E-04 

14 

86400.0 

3.5248 

0.1175E-04 

155.6485 

0.7361E-04 

15 

129600.0 

5.1084 

-0 . 4343E-06 

159.0629 

0 . 6081E-04 

16 

172800.0 

5.9240 

0. 1422E-04 

163.7965 

0 . 7615E-04 

17 

216000.0 

6.7804 

0 . 1 130E-04 

164.9938 

0 . 71 98E-04 

18 

43200.0 

0.3364 

O.OOOOE+OO 

154.4359 

O.OOOOE+OO 

19 

172800.0 

0.5223 

O.OOOOE+OO 

149.6579 

O.OOOOE+OO 

20 

0.0 

7.5385 

0 . 4305E-04 

159.0714 

-0.6106E-05 

21 

172800.0 

9.7906 

0.2576E-04 

157.7525 

0. 1597E-04 

nn 

0.0 

6.4738 

0 . 2285E-04 

158.6602 

-0 . 4225E-05 

23 

43200.0 

7.2620 

0 . 2001E-04 

159.0805 

-0 . 3449E-04 

24 

86400.0 

8.2195 

0 . 3025E-04 

159.2360» 

0.6079E-05 

cr 

- 

172800.0 

11.4950 

0. 1706E-04 

158.0335 

-0 ♦ 8S65E-05 

26 

86400.0 

4.7811 

0 . 3463E-04 

150.2000 

-0.2200E-04 

->7 

216000.0 

13.4423 

0 . 7618E-04 

151.3719 

-0.2750E-04 

28 

43200.0 

6.6095 

-0 . 5367E-04 

155.0578 

0 . 1522E-04 

t’ 

129600.0 

-0.9081 

-0 . 7107E-04 

157.1149 

0 . 2712E-04 

216000.0 

-6.2626 

-0.9251E-04 

159.4644 

0. 2897E-04 

20 
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The  individual  track  files  are  contained  in  both  synthetic  data 
bases.  As  Table  2-6  indicates,  the  only  difference  in  data  bases 
ORT1A  and  ORT1B  is  the  measurement  accuracy  of  the  sensors.  The 
simulated  sensor  accuracy  in  data  base  ORT1A  is  ten  times  greater 
than  that  in  ORT1B. 

The  test  cases  involving  simulated  data  were  run  on  a  DEC  10 
computer. 


2.  2  Problem  Structure  for  Simulated  MSI  Scenarios 

The  basic  structure  of  the  k-track  correlator  is  that  of  an 
integer  program  that  operates  on  a  set  of  potential  correlations  to 
select  the  most  accurate  picture  possible  of  the  surveillance  area. 

The  set  of  potential  correlations  for  Run  #1  selected  by  the 
backtracking  Kalman  filter  process  described  in  Section  1.2.  1  is 
shown  in  Figures  2-2  through  2-12.  Each  figure  indicates  the  track 
file  data  points  involved  in  that  particular  correlation.  Although  the 
parameter  file  setup  allowed  as  many  as  ten  tracks  to  be  correlated 
simultaneously,  the  various  tests  employed  during  the  backtracking 
process  selected  only  2-  and  3-track  correlations. 

Three  basic  types  of  errors  can  occur  in  the  process  of  form 
ing  potential  k-track  correlations: 

Type  1  Error:  the  tracks  included  in  the  correlation  are 
correct  (in  that  they  correspond  to  the  same  ship),  but  not 
all  tracks  from  the  same  ship  were  included. 

Type  2  Error:  the  correlation  includes  tracks  from  more 
than  one  ship. 

Tvpe  3  Error:  failure  to  detect  a  true  correlation. 


Potential  correlation  #1  (ORT1A/PAR 1  A) 


pxHtiW 


Figure  2-3.  Potential  correlation  #2  (ORT1  A/PAR  1  A). 


re  2-4.  Potential  correlation  #3  (ORT1A/PAR  1A) 


Figure  2-6.  Potential  correlation  #5  (ORT1A/PAR  1A). 


Potential  correlation  #6  (ORT1A/PAR 1  A) 
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Figure  2-10.  Potential  correlation  #  9  (ORT1A/PAR  1  A). 


Table  2-7  below  summarizes  the  results  of  the  backtracking  Kalman 
filter  on  a  correlation-by-correlation  basis.  Note  that  both  Type  1  and 
Type  2  errors  are  mixed  in  with  the  set  of  correct  correlations.  No 
Type  3  errors  were  noted  in  either  run. 

Table  2-7.  Summary  of  potential  correlations 
constructed  during  Run  #1. 


Correl 

Number 


Data  Point 


Time 

Number 

.  0 00 00 E +  00 

. OOOOOE+OO 

20 

. 12000E+02 

23 

. 24000E+02 

24 

.  48000E+02 

25 

.  48000E+02 

2:l. 

♦00000E+00 

20 

♦00000E+00 

22 

. 12000E+02 

23 

♦ 24000E+02 

24 

, 43000E+02 

21 

♦48000E+02 

25 

♦00000E+00 

20 

♦  2400 0E +02 

7 

. 36000E+02 

r- 1 

w 

. 48000E+02 

2:i 

. 46000E+02 

w 

♦24000E+02 

la 

. 2 t000E+02 

Comments 

Correct.  Correlates  Tracks 
#7  and  #8. 


Same  as  above.  Note  that  data 
points  were  reordered  due  to 
coincidental  time  of  report. 


[Incorrect  (Type  2  error). 


Incorrect  (Type  1  error).  Cor¬ 
relates  Tracks  #3  and  #5,  but 
leaves  out  Track  #4. 


. 48000E+02 
<  6'.’OOOE-K>2 


*  0  1 0  0  0 E  +  00 

. 24000E+02 

2  a  ( v  02 

#.  •  •  !.  •  ’  I  *  '  *  *  "y  ’  • 


Correct.  Correlates  Tracks 
#3,  #4,  and  #5. 


I 


Table  2-7.  (Continued) 


p- 

Correl 

Data  Point 

1 

r 

Number 

Time 

Number 

Comments 

i 

| 

i 

i 

.  oo  oE+oo 

. 2-000E+02 

10 

1 1 

Incorrect  (Type  1  error).  Cor- 

1 

• 24000E+02 

1  -9 

relates  Tracks  #4  and  #5,  but 

. 36000E+02 

.1.2 

leaves  out  Track  #3. 

♦ 36000E+ 02 

15 

« 48000E+02 

13 

• 48000E+02 

16 

f  » 

, 60000E+02 

17 

>■ 

7 

. 00000E+00 

1 0 

Same  as  Correlation  5.  Note 

. 24000E+02 

l :: 

that  data  points  were  reordered 

.24000E+02 

/ 

due  to  coincidental  time  of 

.24000E+02 

14 

report. 

> 

. 36000E+02 

12 

♦ 36000E+02 

8 

.36000E+02 

15 

* 

. 48000E+02 
. 48000E+02 
.48000E+02 

13 

9 

16 

(  •  » 

♦  60000E+02 

1 7 

\ 

8 

, OOOOOE+OO 

j  0 

Incorrect  (Type  1  error).  Cor- 

I 

. 24000E+02 

:! ;! 

relates  Tracks  #3  and  #4,  but 

| 

. 24000E+02 

leaves  out  Track  #5. 

Sr 

. 36000E+02 

.w 

1  1 

. 36000E+02 
,  -xSOOOE+02 
♦ 48000E+02 

.!.  'S 

I 

9 

♦ 12000E+02 

1 2000 E +02 

Correct.  Correlates  Tracks 

1 

.  24000E+02 

* 

#1  and  #2. 

.36000 E+ 02 

1 

..  48000E+02 

i, 

.  60000E '  + 02 

i' 

•  » 

•/• 

■  ~~ 

Incorrect  (Type  2  error). 

-  3i>000£+02 

v  ovOOOt  ••  0  2 

1 

‘i  vG0:;"  -f  v*2 

Same  as  Correlation  9.  Note 

.  a  k o  ji: £t02 

that  data  points  were  reordered 

<•  lJ  xl  '■ )  0 ' f  f:: '  0  *: 

due  to  coincidental  time  of 

> 

| 

2-CtOOOEtO  2 

-  o i:;  •;■  :• 

r: 

report. 

Note  that  in  Table  2-7  each  individual  correlation  consists  of 
combinations  of  complete  track  files.  For  example,  correlation  10  con¬ 
sists  of  data  points  {l,  26,  2,  27),  and  we  don't  see  correlations  such  as 
{1.  26,  2).  This  particular  problem  structure  is  the  defining  character¬ 
istic  of  the  track-to-track  correlation  problem.  It  arises  because  it  is 
assumed  that  complete  tracks  (e.g.,  {l,  2}  and  {26,  27))  are  constructed 
properly  at  the  sensor  level.  (By  contrast,  in  a  sensor-level  multi¬ 
target  tracking  problem,  the  key  problem  is  to  decide  whether  or  not 
data  points  1  and  2  are  the  same  ship  [1].  ) 

The  structure  of  the  0-1  integer  program  for  k-track  correlation 
can  be  illustrated  by  converting  the  potential  correlations  for  Run  #1  into 
the  form  of  Equations  (11)  through  (13),  repeated  here  for  convenience: 

max  d*? 
subject  to 

Bv  <  1 
v  binary  . 


For  Run  #1,  this  problem  has  the  following  form: 


Track  File  Numbers 


I 


I 


t 


l 


I 


» 


» 


» 


I 


The  vector  d  has  elements  which  correspond  to  the  log- 
likelihood  function  (Equation  (8) )  for  each  potential  correlation. 

Thus  the  first  element  of  d  is  the  numerical  value  of  tnp(u>)  evalu¬ 
ated  at  the  last  data  point  in  Correlation  1. 

Note  that  B  is  a  10x11  matrix  (number  of  track  files  x 
number  of  potential  correlations).  B  is  made  up  of  columns  of  zeros 
and  ones,  each  column  representing  one  particular  correlation.  For 
example.  Column  1  in  B  represents  Correlation  1,  which  consists 
of  Tracks  #7  and  #8  (the  7th  and  8th  elements  of  the  column  are 
set  to  1 ). 


Closer  examination  of  B  reveals  that  the  problem  can  be 
decomposed  into  a  set  of  independent  subproblems.  By  appropri¬ 
ately  permuting  columns  of  B,  we  have  the  equivalent  matrix 


Subproblem  1 


B'  = 


Thus  B'  is  a  decomposition  of  the  original  intersensor  correlation 
problem  into  two  subproblems,  with  one  coupling  constraint. 
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The  structure  of  the  problem  is  further  reduced  by  column 
i 


elimination.  Let  b  be  the  i-th  column  of  B.  Then  the  j-th  column 
of  B  can  be  eliminated  if 


bl  =  b> 
and 


d.  2d., 
1  J 


where  d.  is  the  i-th  element  of  the  cost  vector  d.  Furthermore,  the 
1 


k-th  column  can  be  eliminated  if  d^  *  0.  Finally,  any  row  of  zeros 
can  be  eliminated.  Carrying  out  this  process  for  B1,  we  arrive  at 
an  integer  program  with  maximum  decoupling  and  minimum  size 
(note  that  the  coupling  constraint  was  eliminated,  and  that  Subprob¬ 
lem  1  was  further  decomposed): 


Subproblem  la 


d"  = 


126. 

1 

0 

0 

0 

0 

0 

95. 

1 

0 

0 

0 

0 

0 

187. 

,  B"  = 

0 

1 

1 

1 

0 

96. 

0 

1 

1 

0 

1 

0 

174. 

0 

0 

1 

1 

1 

0 

130. 

0 

0 

0 

0 

• S)  * 

1 

0 

0 

0 

0 

0 

1  - 

Subproblem  lb 


Subproblem  2 


Application  of  these  methods  markedly  decreases  run  time, 
which  depends  exponentially  upon  the  size  of  B. 
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2.3 


Algorithm  Characteristics 

The  truth  models  and  results  of  integer  program  processing 
for  Runs  #1  and  2  are  shown  graphically  in  Figures  2-13  through 
2-16.  As  the  figures  indicate,  in  Run  #1  the  k-track  correlation 
algorithm  performed  correctly.  In  Run  #2,  with  sensors  one  order 
of  magnitude  less  accurate  (ratio  of  standard  deviation  of  noise),  one 
Type  1  error  was  made. 

Note  that  the  integer  program  is  constrained  to  select  from 
among  the  elements  of  potential  correlation  set  D  during  the  process 
of  forming  a  complete  picture  of  the  surveillance  area.  Since  Type  1, 
Type  2,  and  Type  3  errors  can  occur  during  the  formation  of  the  set  D, 
the  integer  program  is  subject  to  errors  of  the  same  type. 

The  total  run  time  requirement  of  the  k-track  correlation 
algorithm  depends  upon  the  run  times  of  the  two  individual  segments 
of  the  code:  (1)  the  construction  of  potential  correlations,  and  (2)  the 
integer  program.  Of  the  two  segments,  the  integer  program  is  the 
most  sensitive  to  program  size.  As  Figure  2-17  indicates,  the 
results  for  ten  test  cases  indicate  a  reasonable  growth  in  the  back¬ 
tracking  Kalman  filter  run  time  as  the  problem  size  increases.  As 
Figure  2-18  indicates,  the  integer  program  exhibits  a  sustained 
exponential  increase  in  run  time  as  problem  size  increases.  The 
underlying  problem  structure  is  such  that  the  decoupling  procedure 
mentioned  in  Section  2.  2  above  is  critical  for  MSI  problems  where 
extremely  large  numbers  of  potential  correlations  are  found. 

In  several  cases,  a  suboptimal  solution  was  returned  by  the 
integer  program  due  to  run  time  constraints  placed  on  the  code.  In 
those  cases  where  the  maximum  time  limit  resulted  in  a  possibly  sub- 
optimal  solution,  good  accuracy  was  obtained.  Thus  it  appears  that 
good  suboptimal  solutions  to  the  integer  program  can  be  obtained  within 
a  reasonable  period  of  time  for  the  types  of  data  bases  discussed  in 
this  paper. 
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All  Correlations 
Correctly  Formed 
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WEST  LONGITUDE 

Figure  2-14.  Correlated  tracks  produced  by  integer  program 
(ORT  1  A/PAR  I  A). 


Figure  2-15.  Actual  tracks  contained  in  synthetic  data  base  ORT1B 


3.0  CONCLUSIONS 


This  paper  summarizes  a  study  which  consisted  of  essentially 
two  parts: 

•  initial  development  of  a  sophisticated  intersensor 
(track-to-track)  correlation  algorithm 

•  preliminary  analysis  of  the  algorithm's  characteristics 
when  applied  to  synthetic  data. 

The  primary  finding  of  the  paper  is  that  an  accurate  picture  of  a  large 
ocean  surveillance  area  can  be  constructed  automatically  by  the  k-track 
correlation  algorithm.  This  finding  is  substantiated  by  the  analysis  of 
synthetic  data  discussed  in  Section  2. 

The  structure  of  the  algorithm  is  that  of  a  Bayesian  decision 
process,  which  by  its  very  nature  may  produce  some  errors  in  the 
analysis  of  any  particular  data  base.  As  noted  above,  these  errors 
were  minimal  for  the  analysis  performed  herein.  There  is  a  critical 
need  for  parameter  tuning  against  truth  models  before  the  algorithm 
can  be  trusted  in  an  operational  situation.  Thus  it  is  clear  that  sub¬ 
stantial  work  remains  before  the  promising  results  obtained  during 
this  short  preliminary  study  can  be  broadly  applied. 

In  the  area  of  track-to-track  correlation  algorithm  develop¬ 
ment,  there  remains  a  number  of  important  issues  yet  to  be  resolved. 
A  few  of  these  are  as  follows: 

•  How  large  a  problem  can  the  integer  program  handle 
in  real  time  after  the  decoupling  procedure  mentioned 
in  Section  2.  2  is  implemented? 

•  Using  the  k-track  correlation  algorithm  as  a  bench¬ 
mark,  how  accurate  are  simpler  suboptimal  correla¬ 
tion  algorithms  (such  as  the  pairwise  algorithm  men¬ 
tioned  in  Section  1. 2)? 


■M—imwnl*»JWI-y- d'i>  •>» 


•  What  is  the  impact  on  algorithm  accuracy  of  Kalman 
filter  mismatch  (i.  e. ,  a  mismatch  between  actual 
and  assumed  sensor  accuracy)? 

•  What  is  the  feasible  operating  regime  for  track-to- 
track  correlation  algorithms  in  terms  of  ship  density? 

•  How  sensitive  is  algorithm  accuracy  to  intersensor 
bias?  How  accurately  can  inter  sensor  alignment  be 
carried  out?  Can  we  "bootstrap"  the  alignment  as 
track-to-track  correlation  is  carried  out? 

•  What  is  the  proper  balance  between  real-time  track- 
to-track  correlation  accuracy  and  the  load  on  sur¬ 
veillance  network  communication  links  (how  much 
data  is  enough  for  a  given  surveillance  area)? 

This  partial  list  of  important  MSI  issues  yet  to  be  resolved 
can  of  course  be  substantially  expanded.  Hopefully,  the  results  of 
this  short  preliminary  study  of  one  specific  correlation  algorithm 
will  resolve  some  of  the  complex  issues  surrounding  automatic, 
accurate,  reel-time  ocean  surveillance. 
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